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CALCULATION OF NONLINEAR CONICAL FLOWS 
BY THE METHOD OF LINES 


By E. B. Klunker, Jerry C. South, Jr., and Ruby M. Davis 
Langley Research Center 

SUMMARY 

A computational technique, called the method of lines, is developed for computing 
the flow field about conical configurations at incidence in a supersonic flow. The method, 
which makes use of the self -similarity property, is developed for the nonlinear flow 
equations. The basic idea is to discretize all but one of the independent variables in the 
partial differential equations so that a coupled system of approximate, simultaneous, 
ordinary differential -difference equations is obtained. Initial values of these differential - 
difference equations are determined from the shock relations after the shock shape is 
estimated or otherwise specified. The system of equations is integrated numerically 
and an iterative process is utilized for adjusting the shock shape to satisfy the boundary 
condition of flow tangency on the body. 

The method has proved to be an efficient and versatile procedure for constructing 
the numerical solutions to conical flow problems. It has been successful in computing 
the flow about circular and elliptic cones at conditions where small regions of supersonic 
cross flow develop and for the conical delta wings where the region of supersonic cross 
flow is extensive. The calculations made for circular and elliptic cones as well as for 
the compression side of various conical delta wings are in good agreement with experi- 
ment except in regions where viscous effects become important. 

INTRODUCTION 

In 1935 Busemann (ref. 1) introduced the concept of a general conical flow field as 
one in which the fluid properties are constant along any ray emanating from a common 
point in the flow. Solutions for such self -similar conical flows are of great importance 
to the aerodynamicist since (1) significant regions of the flow about many practical con- 
figurations are conical, or nearly so; (2) conical bodies and wings are the simplest class 
of three-dimensional shapes and thereby provide "benchmark" cases for both theoretical 
and experimental studies in supersonic and hypersonic flow. 

Although the self-similarity property of conical flow allows the reduction of the 
problem from three to two space dimensions, the analyst finds himself confronted with a 
formidable free -boundary problem for nonlinear partial differential equations of elliptic 



or mixed type. Hence until the last decade, most conical solutions have been obtained 
only for the simplest cases or after linearization or other approximations to the equa- 
tions. However, recent advances in speed and storage of digital computers have spurred 
the development of numerical solutions of the full nonlinear equations, to the point where 
solutions to very general conical flow problems can be obtained in a few minutes. 

A particularly efficient numerical technique for solving conical flow problems has 
been reported in references 2 to 4. The method is semidiscrete, wherein one independent 
variable is discretized while the other remains as a continuous variable. This method is 
referred to as the method of lines, to distinguish it from grid or network computations 
where all independent variables are discretized. 

The method of lines is "direct” in the sense that the body shape is given and is one 
of the bounding coordinate surfaces; yet, the shock wave is another bounding coordinate 
surface, and the governing differential equations are solved by integrating inward from 
the shock. Thus, in that respect, the method is like the inverse methods. The technique 
employed for solving this free boundary problem has three distinguishing features: 

(1) the coordinate transformation which maps the region between the shock and body onto 
a rectangle, (2) the solution of the equations by the semidiscrete method of lines, and 
(3) an iteration procedure for satisfying the boundary conditions. None of these features 
are new; yet when combined, they prove to be an efficient means of solving free boundary 
problems such as the supersonic blunt -body problem or conical flows. The basic idea of 
the method of lines is to discretize all but one of the independent variables in the partial 
differential equations so that a system of approximate, simultaneous, ordinary, 
differential -difference equations is obtained. Initial values for the system of equations 
are estimated, or otherwise specified, and the system of equations is integrated numer- 
ically. An iterative process is utilized to satisfy the boundary conditions; thus, the initial 
values are subsequently altered and the equations are again integrated. The success of 
the method of lines as a computational tool hinges upon (1) formulation of the problem in 
a form which requires relatively few lines, (2) use of an efficient integration routine that 
yields good accuracy with relatively large integration steps, and (3) development of an 
efficient interative process to satisfy the boundary conditions. The first requirement is 
largely met through the choice of the coordinate system and the second can be satisfied 
with any of a number of integration schemes such as a fourth-order Runge-Kutta method. 
The computational time and the utility of the method depend to a large part on the itera- 
tive process. 

The present paper expands upon the material in reference 4; further details and 
refinements of the method are presented, together with numerous applications to a 
variety of conical flow problems. Comparisons of the present calculations with other 
theories and experiment are given for circular and elliptic cones, and conical delta wings. 
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The stream velocity vector lies in a plane of symmetry for all the configurations; how- 
ever, this restriction is not a limitation of the method itself. Some of the more theoret- 
ical aspects of conical flow are touched on, such as the inviscid entropy layer with the 
attendant steep gradients adjacent to the surface, and the nodal-type singularities of the 
cross -flow streamline patterns. 

To aid the reader, numerous headings and subsections are employed. A separate 
section "Background" is included which cites most of the recent work in nonlinear coni- 
cal flow theory, in particular, the related work in the U.S.S.R. which seems to have gone 
relatively unnoticed. 


BACKGROUND 

In this section, a review of past work in nonlinear supersonic conical flow theory 
is given. No effort has been made to consider the large body of literature which con- 
cerns linearized conical flow theory. 

Nonlinear Conical Methods 

The earliest treatment of nonaxisymmetric conical flow was given by Stone (ref. 5) 
together with the numerical computations carried out under the direction of Kopal (ref. 6), 
where the flow about circular cones at small incidence was constructed as a perturbation 
about the axisymmetric nonlinear Taylor -Maccoll solution (ref. 7). Ferri (ref. 8) rec- 
ognized the singularities of these conical flows, which were not accounted for in the 
Stone solution, and discussed the general features of the streamlines. The many analyt- 
ical papers published since (refs. 9 to 11 and many other papers referenced in these 
works) have been concerned largely with the construction of solutions to conical flows by 
means of matched asymptotic expansions. These papers have concentrated primarily on 
the theoretical development; consequently, there has been relatively little computational 
work presented. 

Two basic approaches are available for the numerical development of exact non- 
linear conical solutions 

(1) Distance -asymptotic methods, where some initial distribution of the flow vari- 
ables and shock-wave shape is used near the apex as initial values for continuing the 
calculation downstream by some three-dimensional computation scheme. The calculation 
proceeds until conical similarity conditions are sufficiently satisfied. 

^Exact in the sense that the only approximation made is the reduction of the gov- 
erning partial differential equations to ordinary differential equations or algebraic equa- 
tions by using finite -difference expressions for the derivatives with respect to one or 
more of the coordinates. 
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(2) Methods which invoke the conical self -similarity and thereby reduce to two the 
number of independent variables are referred to simply as "conical" methods for brevity. 

Both general approaches have their merits. The distance -asymptotic techniques 
develop the solution as a well -posed initial-boundary problem for equations of hyperbolic 
type, and convergence is "almost" guaranteed from both physical and theoretical consid- 
erations. However, to achieve a satisfactory solution in many problems where a fine 
mesh is needed, these methods require a large amount of computer storage and time. 

The conical methods reduce the problem to one in two dimensions but in the more diffi- 
cult form of a free boundary problem for equations of elliptic or mixed type. In fact, 
many of the conical methods are similar to methods used for solving the blunt-body 
problem. 

The methods of references 12 to 20 are examples of the distance -asymptotic 
method. References 12, 13, and 14 considered circular cones at angle of attack, refer- 
ences 15 and 16 included cones of elliptic cross section, and reference 17 presents cal- 
culations for the compression side of conical delta wings with the shock wave attached 
not only at the apex, but also along the swept leading edges. More recently, a three- 
dimensional characteristic method has been developed to compute the flow about some 
delta -wing configurations, also with an attached leading-edge shock (ref. 18). In refer- 
ences 19 and 20, a method is presented which is essentially a distance -asymptotic 
method, and which treats shock waves not as discontinuities, but as rapid but continuous 
compressions. This method appears to be most promising for problems in which com- 
plicated embedded shock patterns may occur. Reference 19 presents results not only 
for circular cones at incidence, but also for a conical wing -body configuration. Refer- 
ence 20 includes results for both the expansion and compression sides of planar delta 
wings with shock attached at the leading edge. 

Examples of conical methods appear as early as 1929, when Busemann (ref. 21) 
constructed the axisymmetric conical flow by numerical -graphical construction in the 
hodograph plane. Reference 22 (pp. 526-536) cites many of the approximate and exact 
conical methods documented up to about 1964; therefore, they are not all discussed in the 
present paper. Most recent conical methods are of the "inverse" type, in which a simple 
analytic function is assumed for the conical shock wave, and the governing partial differ- 
ential equations are solved by marching inward until some body shape is obtained. Var- 
ious inverse methods are reported in references 23 to 27. These methods have not been 
successful in constructing solutions for body shapes which produce a shock wave 
requiring many parameters for an adequate description; only circular -cross -section 
cones at incidence have been amenable whereas solutions for elliptic cones have been 
obtained only painstakingly. Other conical methods have used the method of relaxation 
in regions where the cross flow (velocity component normal to a conical ray) is subsonic, 
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and the governing differential equations are of the elliptic type; a two-dimensional method 
of characteristics is used in supersonic cross -flow regions where the equations are 
hyperbolic. Such approaches are developed in references 28 to 30 in which flat delta 
wings with attached leading-edge shocks are considered. It is not known whether these 
methods can be easily coded for efficient machine calculation for nonflat conical wings. 

Method of Lines 

The method of lines is somewhat similar to the better -known (in the United States) 
method of integral relations (ref. 31) in that both methods are semidiscrete. The method 
of integral relations, in its usual formulation (refs. 31 to 35), requires new algebraic 
development for each higher approximation (that is, more lines) and the equations grow 
more complex; whereas the method of lines system is written recursively with arbitrary 
number and spacing of lines. 

Liskovets (ref. 36) presents a general review of the method of lines as applied to 
linear equations of elliptic, parabolic, and hyperbolic types. It is of interest to note that 
the method apparently dates from 1939 with an application to the solution of the Laplace 
equation by M. G. Slobodyanskii. An extensive bibliography (covering the work done in 
the U.S.S.R. up to 1965) is presented in reference 36. 

In reference 37, the work of Telenin and his coworkers is cited, in which they used 
the method of lines for constructing numerical solutions to the axisymmetric, supersonic 
blunt-body problem. That work was extended to the three-dimensional blunt -body prob- 
lem and reported in reference 38.^ More details of the preceding work are given in 
reference 39. 

Makhin and Syagayev (ref. 40) recognized the difficulties associated with Syagayev's 
earlier inverse method (ref. 27) and applied the transformation discussed previously in 
which the body and shock become bounding coordinate surfaces. Their procedure was 
essentially the method of lines. They used a first-order Euler method to integrate 
toward the body and to maintain accuracy, a small integration step size was required 
(l/64th of the local shock-layer thickness). The modification allowed them to solve the 
direct elliptic-cone problem, which was intractable by the earlier inverse approach. 

Bazzhin and Chelysheva (ref. 41) applied the method to conical bodies at large 
angles of attack, where regions of supersonic cross flow always occur. Their procedure 
was very similar to approaches to the blunt-body problem; they used the method of lines 
on the windward, high-pressure side of the conical body up to, and beyond, the region 
where the cross flow becomes supersonic. A conical, two-dimensional method of char- 
acteristics was used to continue the solutions of the supersonic cross-flow region. More 

2The word "Attached" in the title of reference 38 is an obvious error in transla- 
tion, and should read "Detached.” 
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results, with emphasis on the characteristics solutions, are given in reference 42. An 
interesting result in reference 42 was that in some instances, it was possible to continue 
the characteristics solution through the leeward plane of symmetry of elliptic cones at 
large incidence; thus, symmetry conditions were violated and indicated the impossibility 
of flow without embedded shocks in the general case. Neither reference 41 nor 42 sug- 
gests using the method of lines for the small-to-moderate incidences in which the cross 
flow is everywhere subsonic, nor was any mention given to the application in 
reference 40. 

Jones (ref. 2) reported a method similar in many ways to the present procedure. 

He obtained solutions for circular and elliptic cones and a conical body with a four- 
parameter, smooth cross-section contour having both concave and convex portions. The 
method was shown to be accurate and efficient, and he was able to experiment numerically 
with some of the more theoretical questions concerning conical flows such as the "lift-off” 
of the vortical singularity (refs. 8 and 9). Ndefo (ref. 3) has applied the Telenin approach 
to the computation of conical flows. His development of the method of lines, which is 
similar in many respects to that of Jones and the present work, has been applied only to 
the computation of the flow about circular cones. In reference 4, the method of lines was 
applied to circular and elliptic cones, and to conical delta wings with attached leading- 
edge shocks. The present paper develops the method further and presents additional 
computations. 


SYMBOLS 


a semimajor axis of ellipse 

b semiminor axis of ellipse 

C D drag coefficient 

C-^ lift coefficient 

C pitching-moment coefficient about X-axis 

Xil jX 

Cp pressure coefficient 

Cy normal-force coefficient 

C z axial -force coefficient 
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c 


ratio of speed of sound to critical speed. 



e T ,e^,e r unit vectors along r, r j, and r 
h scale factor for r coordinate at r = 1 


U 

U,k 

K 



indices indicating line number 

unit vectors along X-, Y-, and Z -directions 

mean curvature of surface rj = Constant 

mean curvature of conical body given by equation (3) 


free -stream Mach number 


cross -flow Mach number, 



number of lines 


p pressure, referenced to product of stream density and square of critical 

speed 

r distance along ray 

S entropy 

u,v,w components of velocity in r-, rj-, and r-directions, respectively, referenced 

to critical speed 

u c ,v c’ w c cylindrical components of velocity in axial, radial, and azimuthal directions, 
respectively, referenced to critical speed 

V total velocity, referenced to critical speed 
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N|Kj NIX 


X,Y,Z 


a 

13 

y 

e 



V 


V s 

9 

e o 

A 

l 

P 

cr 


Cartesian coordinates 


nondimensional spanwise coordinate, 

x max 


angle of attack or incidence 

shock angle 

ratio of specific heats 

convergence criterion on normal velocity at conical surface 


angle measured in plane normal to body from ray on surface of body to 
ray in field 

value of 7] at shock 

conical apex angle of body in horizontal plane Y = 0 

conical apex angle of body in plane of symmetry containing free-stream 
velocity vector 

sweep angle 

transformed arc length variable 

density, referenced to stream density 

angle between shock normal and 77 -direction (fig. 38) 

arc length along intersection of unit sphere with conical body 
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polar angle in cylindrical coordinates (fig. 2) 


xjy yaw angle 

Subscripts: 

n normal to leading edge 

o body surface coordinate 

r,t],T indicates directions associated with these components 

00 free stream 


METHOD 

Conical Coordinates 

The equations governing supersonic, inviscid flow of an ideal gas are written in a 
body -oriented, orthogonal, conical coordinate system (r,??, t) as developed in reference 43 
where r is the distance along a conical ray, r] is the angle measured from the body 
surface to the ray in a plane r = Constant, and r is a measure of the arc length^ along 
the intersection of the body surface with a sphere of radius r centered at the body apex. 
(See fig. 1.) The body is the conical surface 77 = 0, and the contour along which t is 
measured is a plane curve only in the special cases of a circular cone or a flat delta 
wing. This coordinate system has the advantage that the associated velocity components 
in the r-, rj-, and r-directions are the natural conical components u, v, and w; the 
v-component is zero at the surface 77 = 0 and yv% + w^, the magnitude of the velocity 
component normal to a conical ray (hereafter called the "cross-flow" component), gov- 
erns the type (elliptic or hyperbolic) of the partial differential equations of the conical 
flow. F urtherm ore, certain singularities appear where the cross-flow component van- 
ishes = 0), and it is extremely important that such points be recognized and 

interpreted correctly. These coordinates prove to be advantageous from the computa- 
tional point of view as well except for bodies with concave curvature. 


3ln reference 4 both £ and r are arc length variables where the duplicity arises 
because of a coordinate transformation. In the present development the variable r is 
the arc length and = £(r) is a transformed coordinate which can be selected to vary 
the line spacing. 
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Figure 1.- Body-oriented, conical coordinates. 

Differential Equations 


With the conical similarity = Oj, the partial differential equations which 

describe the flow involve two independent variables, 77 and r. Consequently, the solu- 
tion can be developed on the spherical surface r = 1. The equations of continuity, 
momentum, and conservation of entropy S along a stream surface in these coordinates 
(ref. 43) become 

i — + h — + — + 2 uh - hKv = 0 1 

P Dt dV 9r 

— - h(v 2 + w 2 ) = 0 
Dt 

^ + huv + hKw 2 = 0 > (1) 

Dt P 97? 

■5^ + ~ + huw - hKvw = 0 

Dt P 9 t 

ds = o 

Dt J 

where the total derivative is 
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_D 

Dt 



The scale factor h for the r coordinate and the mean curvature K of the surface 
r] = Constant at r = 1 are 


h = cos t? - Kg sin r] 

K = -5^4( sin ’ ,tK B cos,) ) 


where Kg is the mean body curvature. 

Transformation to a rectangular region. - The integration of the system of equa- 
tions, as discussed previously, is facilitated by a coordinate transformation which maps 
the region bounded by the shock and the body into a rectangular domain as shown in 
figure 2. For this purpose the transformed variables are taken as 



£ = !(t) 


where 77 = t](t) is the shock surface. Thus, £ = 0 on the body and ? = 1 on the 
s 

shock. The transformation to the new independent variable £ permits a local 



Physical (t, n ) plane Transformed (E;, c; ) plane 

Figure 2.- Layout of computational lines. 
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stretching of the arc length coordinate. The chain rule gives 


_a_ = j__a_ 

3 V V s a? 

±-l (± 

9 t ti s d£ at) 


where 




= d| 

dT 


For a perfect gas the derivative of the density in the continuity equation is replaced 
by one involving the pressure through the use of the equation for the conservation of 
entropy along a stream line. Thus, 


Dp _ _1_ Dg 
Dt c 2 Dt 


where c is the ratio of the speed of sound to the critical speed. Then the equations (1) 
in the new coordinates can be written as 


If = §( F 1* - F 2 h + F 3 C£ 


^s\ 
T d *) 




av _ l/h ap 

a? "s\p aC 



aw _ l L, d?? s\ 1 ag 
a? g \* T d| /P a? 



au _ f 

a? e 


( 2 ) 


as _ * w as 

a ? ~ n s^r g 9 £ 


where 


g = hv - w££ 


5_s 

T d£ 
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D = h 2 


« T 




F l = ’'43l! t£ r|f +2 “ h - hKv ) 


F 2 = ^s 1 


(w£, T + huv + hKw 2 j 

F 3 = ( 7 ? |f + w *t ff + huw - M*™) 

f = ; ?s [h(v 2 + w 2 ) -w| T ||] 

and the entropy is related to the pressure and density by 
S = log^yM^c* 2 £pj 


where 


c*=. 2^1* 2 

V y + 1 (y + i)m m 2 

is the ratio of the critical speed to the free -stream speed. The Bernoulli equation 
— ^ + U 2 + V 2 + w 2 = - - \ 

y - IP y - 1 

provides an additional relation between the flow variables. The velocities are nondimen- 
sionalized with the critical speed, the density with the stream density, and the pressure 
with the product of the stream density and the square of the critical speed. 

The right-hand members of equations (2) involve the unknown shock shape, the 
velocity components, pressure, density, and the derivatives of these quantities with 
respect to the variable 4. The density can be computed from the value of the entropy 
and the pressure in the integration of the system of equations (2) or through the use of the 
Bernoulli equation. 

The problem in the rj,T coordinates is thus transformed to the solution of a system 
of equations in a rectangular domain with coordinates £,|. The shock wave and the body 
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surface are mapped onto the lines ? = 1 and £ = 0, respectively. The origin of the arc 
length t as well as the transformed variable ij are taken in the windward plane of 
symmetry; the terminal value of £ is taken in the leeward plane of symmetry for circu- 
lar and elliptic cones and at the leading edge of the conical wing configurations. 

Geometric Parameters 

The body geometry enters the system of equations (2) through the body curvature 
and the arc length r along the intersection of the conical body with the spherical sur- 
face r = 1. These quantities are evaluated from the equation defining the shape of the 
body. The body is defined in the Cartesian coordinate system shown in figure 1. The 
body axis is in the Z -direction and the angle of incidence a is measured relative to 
the Z-axis. From the conical similarity, the flow is constant along each ray defined 

X Y 

by — = Constant, — = Constant. Let 

y Z Z 



y 


Y 

Z 


and let the subscript o denote values on the body. Then the conical body can be defined 
by an equation of the form 



0 


The mean body curvature at r = 1, which is required for the evaluation of the scale 
factor h and the coordinate curvature K (ref. 44), is 


K 


B 




Gy Gjqj + G x Gyy - 2G x GyG X y 


\/M 3/2 

\ A 2 / 


( 3 ) 


where the single and double subscripts on G denote the first and second derivatives 
of G with respect to the indicated variables x 0 and y 0 and 


A i = 1 + x o 2 + y 0 2 

A 2 - Gx + Gy + ^ x qG x + ^oGy^ 
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The differentials of the arc length at r = 1 and the transformed coordinate £ are 
(appendix A) 


dr = 1 A 2 1 '^ 2 G y " 1 dx Q = A 1 ” 1 A 2 1 ^ 2 G x _ 1 dy 0 

d£ = | T dr 


(4) 


The origin of the arc length is taken in the plane of symmetry; thus, the initial values are 
x Q (0) =0, y Q = y o (0), and r = 0. Integration of these equations then determines the arc 
length in terms of the body coordinates. 

In practice equations (4) are integrated numerically to obtain values of r and £ 
at the leeward plane of symmetry for the circular and elliptic cones or values of r 
and £ at the leading edge of wing configurations. The line spacing is specified in terms 
of £ and equations (4) are again integrated to compute the values of x Q and y Q at 
each line. The quantities ^A^ and ^Ag (appendix A) are the normalizing factors for 
the vectors normal to the spherical surface r = 1 and the conical surface G(x Q ,y 0 ^ = 0 , 
respectively. 

The relation between unit vectors e^, e T , and § r in the 77 -, r-, and r-directions 

and the Cartesian unit vectors i, j, and k is required for the development of the 
shock relations and for the computation of the Cartesian coordinates of points in the flow 
field. The direction cosines are given in appendix A. 


The Method of Lines 

The region of interest in the ?,£ -plane is divided by N lines parallel to the 
£-axis (N + 1 lines for the delta wings with attached leading-edge shocks); the line 
i = 1 is taken in the windward plane of symmetry. The stream velocity vector lies in a 
plane of symmetry for all the results presented; however, this is not a limitation on the 
method itself. It is not necessary that the distance between lines be of equal width since 
the spacing of the lines in the physical 77 , t- plane can be adjusted through the selection of 
the coordinate transformation £ = £(t). Figure 2 illustrates the division in the physical 
77 , r coordinates and the transformed %,£ coordinates with £ = t for an elliptic cone 
with nine lines. Also shown is the cylindrical polar angle 0 ; 0 is -it / 2 in the wind- 
ward symmetry plane and tt/2 in the leeward symmetry plane for conical configurations 
at incidence. For the elliptic cone at yaw, 0 = 0 in the windward symmetry plane and 
0 = 7 r in the leeward symmetry plane. 

At each strip boundary or line, the system of equations (2) is reduced to a set of 
ordinary differential-difference equations by replacing the derivative d/ d£ by finite 
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differences. The derivative of the Lagrange interpolation polynomial (ref. 45, for 
example) was used in the computer program with an equal number of lines on either side 
of the line at which 8/ 8£ is computed; therefore, central differencing is obtained when 
the line spacing is equal. In forming £ -derivatives at lines near the leading edge of 
wings with an attached leading-edge shock (the leading edge is line i = N + 1), the num- 
ber of points used in the derivative formulas is reduced, if necessary, to retain an equal 
number of lines on either side of the line at which 8/a£ is computed. A five -point 
formula was used for most computations. 

The system 'of equations (2) is integrated simultaneously along each line 
i = 1, . . ., N with the derivatives in the right-hand members evaluated with the Lagrange 
formula for derivatives. This polynomial approximation for 8/8£ causes the differen- 
tial equations along any line to be coupled to those along the other lines. The system of 
equations (2) thus becomes a system of 5N ordinary differential equations which are inte- 
grated by a Runge-Kutta or similar method. A fourth-order Runge-Kutta integration has 
been employed in the computer program, and the integration step size was generally in 
increments of -0.1 from the shock to a value of £ of 0.1 and by increments of -0.05 
and -0.025 thereafter, except in calculations where details of the entropy layer are 
sought. Note that the system of equations (2) is of order 4N if the density is computed 
with the Bernoulli equation and the equation for the conservation of entropy is deleted. 

The system of 5N equations has been employed in the computations presented herein. 

Symmetry and Boundary Conditions 

Symmetry conditions.- The Y,Z -plane is the plane of symmetry (fig. 1) which con- 
tains the stream velocity vector for all configurations other than the elliptic cone at yaw. 
The stream velocity vector lies in the plane of the major axis (X,Z -plane) for the yawed 
cone. The origin of the arc length r and the transformed variable £ is taken in the 
windward plane of symmetry. The symmetry conditions require that all the flow quan- 
tities other than w are symmetric and the circumferential component of velocity w is 
antisymmetric about the plane of symmetry. Symmetry is accounted for in the formula 
for the £ derivatives by properly reflecting points (lines) about the symmetry plane. 

For the elliptic cone both the windward (i = 1) and leeward (i = N) lines correspond to 
symmetry planes. For the compression side of delta wings, the line i = 1 is a sym- 
metry plane. 

Flow tangency at surface.- The natural coordinate system has been employed in 
part to simplify the form of the boundary conditions on the body. The condition of flow 
tangency on the body surface £ = 0 requires that the normal component of velocity v 
must vanish. Thus, 

v(0,t) = 0 (5) 
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Since the system of equations (2) is integrated along discrete lines, the boundary condi- 
tion (5) is satisfied only on these lines. 

Shock -wave conditions. - The shock -wave shape is initially unknown and must be 
determined through an iteration procedure as described subsequently. The flow variables 
behind the shock are found from the shock jump conditions. (See appendix B.) 

Attached shock at wing l eading edge . - In the present paper the method is applied to 
delta -wing configurations with the shock wave attached along the sharp leading edges. 
Under these conditions, the slope of the shock wave and the flow variables at the leading 
edge are determined directly from the shock conditions and the wing geometry. (See 
appendix B.) Hence, the flow at the leading edge is completely determined without 
recourse to the integration of the system of differential equations. The leading-edge 
condition takes the place of the symmetry properties for closed bodies in the leeward 
plane of symmetry. 


Determination of the Shock Shape 

The form of the shock -wave cross section is given by the unknown function 

t] = V S ( T )- K V s and dv s /d T were known, all the information concerning the shock- 

wave geometry would be known, and the values of the functions p, p, u, v, and w 

at the shock wave (£ = 1) could be evaluated from the shock jump conditions. These 

values could be used to start the numerical integration at £ = 1 and proceed down to 

the body surface at £ = 0. Only the correct shock function n will cause the flow- 

s 

tangency condition (5) to be satisfied; thus, there must be a relation between the func- 
tion V S ( T ) and the normal component of velocity at the surface, v(0,r). This is the 
basis, then, for determining the correct shock shape. 


Newton iteration for shock shape.- The number of unknown values of n is equal 
to the number of normal components v^(0) = v(0,r^. A Newton-type iteration procedure 
is used for adjusting the N values of ? 7 g . to achieve m.ax|vj(0)| ^ e where e is the 

prescribed accuracy criterion. The steps in the procedure, which are straightforward 
and easily automated, are as follows: 


(1) Assume an initial set of values rj , (i = 1 . . N) based on experience, 

S,1 

approximate solutions, or previously computed cases with conditions close to those 
desired. 


(2) Compute 
derivatives. 


drj 


at each line with the use of the polynomial expression for the 
d?7s i 

(3) With rj ■ and — — -f- from steps (1) and (2), solve the shock jump conditions 
“j 1 d | 


for 


Pi’ Pi’ 


u i’ v i’ 


and wj (i = 1, . . ., N). 
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(4) Use the results of step (3) for initial values to start a numerical integration of 
the system of 5N equations from £ = 1 to £ = 0, and evaluate the surface normal com- 
ponents v^(O). 


(5) Test maxjvj(0)| . 


If 



the shock shape is satisfactory and the 


problem is considered to be solved. Otherwise, the following perturbation cycle is used. 


(6) Perturb each parameter r? g j independently and in order. That is, change 

T] . to (1 + 6)r] ., where 5 is a small number (for example, 10-6), and repeat 
SjJ S ?J 


steps (2) to (4). This procedure results in small changes in each of the v^(0) due to 


the perturbation in r] 


s ,r 


coefficients or partial derivatives 


Hence, the jth column of an 

avj (0) 


drj 


(]' = 1, 


N by N matrix of influence 
. ., N) is generated with each 


s,J 


perturbation. 


(7) Solve the usual first-order linear system 


N 9v.(0) 
J=1 


A7 ?s,j = - v i (0) 


(i = 1, • • N) 


( 6 ) 


to obtain the increments An . required to correct the shock shape and drive all 

s,i 

Vi (0) - 0. 

(8) Use the new shock parameters n . + An . to start a new cycle at step (2). 

S y 1 S J 1 

Note that a complete cycle requires N h 1 integrations: one "pivotal" and 

N perturbations. 

Modified Newton iteration procedure.- The regular Newton iteration procedure 
calls for a reevaluation of the Jacobian influence matrix at the beginning of each new 
cycle (step 6). Since most of the computation time is taken up in this step, a modifica- 
tion which bypasses this step on subsequent cycles (if they are required) was introduced. 
Once the surface velocities are within some preassigned magnitude (considerably larger 
than the final convergence criterion), the components of the Jacobian are not recomputed 
after each pivotal integration; only the right-hand member of equation (6) is updated at 
each cycle. Although more cycles are necessary for convergence, each cycle with the 
constant Jacobian consists of only one integration, and a considerable saving in the total 
number of integrations, and hence computer time, is obtained. 

It should be noted here that the modified Newton procedure can diverge in certain 
cases where the regular Newton procedure converges. The situation is easily visualized 
in the one-dimensional case where the first guess is too far away from the solution and 


18 



the second derivative is large. Figure 3 is an illustration of two sequences of modified 
Newton iterations to determine the zero of a function f: the unprimed sequence starts 


f(x) 



Figure 3.- Converging and diverging sequences of modified Newton iterations. 

at x Q and diverges, but the primed sequence converges. If the slope had been recom- 
puted one more time at x^, the unprimed sequence would also converge. For this 
reason the option is retained to use the regular Newton for more than one cycle if the 
surface velocities are still too large after the first iteration. In nearly all cases com- 
puted by the authors, the modified Newton procedure converged if the max|v^(0)| was in 

the range of 0.03 to 0.05. Examples can be given where the regular Newton procedure 
also diverges. The sufficient conditions for convergence of both procedures are the 
same (ref. 46), but the graphical example of figure 3 demonstrates that the necessary 
conditions (for convergence) are clearly different. 

The shock-wave determination procedure of reference 2 differs somewhat from the 
present one, in that the shock function is given by a trigonometric polynomial, and the 
coefficients are chosen to minimize the sum of squares of the normal components at the 
surface. 

Approximate starting shock shapes.- Usually a very good estimate of the shock 
shape is required for a successful calculation and convergence. The exceptions are cir- 
cular cones at moderate relative incidences, ~~ ~ 0.5, and delta wings at large super- 

° o 

sonic Mach numbers, > 3.0. For most other cases, however, considerable care 
must be exercised in choosing the initial shock shape to start the iterations previously 
described. A poor initial estimate can result in any of several program failures, such 
as negative pressures or vanishing denominators. The latter difficulty is caused by 
excessive supersonic cross flow and introduces characteristic -type singularities in the 
differential equations (that is, when a line £ = Constant is tangent to a conical 
characteristic) . 
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In order to obtain a good initial estimate of the shock shape for such cases as 
elliptic cones at incidence, a "simpler" case is computed first and then the input param- 
eters are changed in a series of steps toward the desired configuration, a new converged 
shock shape being obtained with each change of the input parameters (for example, y, 

M m , 9 q , a, and a/b). The history of the set of r? g . as a function of the set of input 
parameters is incorporated in an extrapolation procedure to predict the new shock shape 
for the new set of input parameters. The converged values of rj . for the initial value 

S,1 

of the parameter are used as the input values corresponding to a small variation of the 
input parameter; for example, the parameter may be changed to 1.05 times its initial 
value. Once the converged solution for this value of the parameter is obtained, the 
inputs T) ■ for a new value of the parameter are computed by linear extrapolation of the 
two previous sets. After three sets of converged 77 1 have been obtained corresponding 
to three values of the input parameter, including the set with the small variation of the 
parameter, a quadratic extrapolation is employed. Such a procedure was found to be 
essential for efficient computation. The procedure is completely automated in the com- 
puter program for circular and elliptic cones, and variation of any of five different input 
parameters is allowed. The program for the conical delta wings has not been so 
automated. 


If, during a sequence of calculations in which a parameter is incremented, the cal- 
culation encounters difficulties (such as negative pressure), the program halves the 
increment and attempts a restart from the last parameter value for which the solution 
converged. If difficulties are encountered again, the parameter increment is halved 
again and this process is automatically repeated until the increment is smaller in mag- 
nitude than a preassigned value, or until successful convergence is achieved. 


The initial estimate for the delta-wing shock shape can be made directly for any 
angle of attack a (up to that causing detachment from the leading edge), sweep angle A, 
and M m . The estimate used is an even function in | which requires only an estimate 
of ri a ■, (the value of rj in the plane of symmetry). The function is contrived to give 
n N+l = ® (tt> e conc fcti° n attachment at the leading edge) and the correct value for 

/d7? A 

— .£ , which is calculated from the shock conditions. 

V d £ /N+l 


The function is 


s ,i 



V 


«i 


5 N+l, 


»N+1 


's,l 


fag I 

\ d ? / 


N+l 


(7) 


where 

tions. 


d 7 s h 

— 5 = - — tan cr 
d| £t 

(See appendix B.) 


and a is evaluated from the wing geometry and stream condi - 

The value of n . used in equation (7) is a tangent-cone 

s , 1 
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approximation, increased by a factor 1.2 to avoid the Mach wave condition for very thin 
wings at small a as follows: 

where (ref. 47) 

AVA oo 

In a few cases, the value of rj 1 was so far off that the required corrections 

s, i 

A ?7 . were sizable, too much "roughness" in the shock shape resulted, and subsequent 

s,i 

iterations failed. It happens that in this event, the first correction for n n is very 

Sj J- 

good so that using that value in equation (7) and restarting the iteration has always been 
successful. 


Extrapolation to Surface 

At the surface, ? = 0, the derivatives du/d? and dw/d? are infinite because of 
the well-known vortical layer adjacent to the surface in the conical flows. (See ref. 8.) 
The derivatives dp/d? and dv/d? however are finite, and this fact allows extrapolation 
of the functions p and v to the surface from ? > 0. In fact, it is only Vj that is 
extrapolated during the integrations in order to evaluate the magnitudes of v^(0) 

(i = 1, . . ., N). When the convergence criterion is met, the pressure p is also extrap- 
olated to ? = 0. 

Corrected isentropic surface values.- The surface entropy is a constant on the sur- 
face, and if the value is known, the isentropic surface density can be calculated as a func- 
tion of the extrapolated pressure and the surface entropy. The Bernoulli equation 
relates p, p, u, v, and w; since v = 0 at the surface and p and p are known, 
there is only one equation with two unknowns, u and w. The other equation needed for 
determining u and w is the differential equation for the u-momentum. This equation 
on the body becomes 


w(0,£) = b T 


9u(0,£) 

db 


(8) 


Substitution of equation (8) into the Bernoulli equation gives the set of nonlinear differen- 
tial equations which are to be satisfied at each line: 


-221. + u. 2 

y - l Pi 1 


+ 



_ y + 1 
y - l 


(i = 1, . . ., N) 


( 9 ) 
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The Lagrange derivative formula is used for evaluating the i; -derivatives, and the set of 
equations (9) is solved by Newton iteration. When convergence is achieved, equation (8) 
gives Wj; hence, the isentropic surface values are completely determined. This pro- 
cedure is essentially the same as that of reference 2. 

The derivatives av/a? and ap/a? at ? = 0 can be evaluated from the corrected 
surface values. The reduction of the continuity equation and the v and w momentum 
equations at £ = 0 gives the relations 


av 

a? 



«2 + „ 

-*> / T 



if - -” S K B PW 


2 


The relations show, as stated previously, that 8v/a£ and ap/a? are finite at the sur- 
face. The finiteness of these quantities is one of the principal merits of the coordinate 
system employed. 

Computation of the surface entropy.- For the circular cone, the surface entropy is 
the value that occurs in the windward symmetry plane, line i = 1. For the elliptic cone 
with the free -stream velocity vector lying in the plane of the minor axis, the surface 
entropy is assumed to be the maximum value at the shock wave. When the free-stream 
vector is in the plane of the major axis, the surface entropy is piecewise constant for 
some flow conditions; that is, it has the value of the windward symmetry plane on the 
surface segment where w(o,r) > 0; and it has the value of the leeward symmetry plane 
on the rest of the surface. The surface entropy for the delta wings with convex surfaces 
is the same as the leading -edge value. 


Stability and Error Growth 

In practical numerical computations by the present method, it is found that there 
exists a maximum number of lines (or minimum spacing between lines) beyond which 
instabilities swamp the solution. This result is by no means surprising, in light of 
Hadamard’s (ref. 48) famous example of the inherent instability of the Cauchy (initial- 
value) problem for Laplace's equation. The example demonstrated that certain types of 
initial data, described by vanishing amplitude and increasing frequency, are magnified 
exponentially as the distance from the initial line increases. No matter how close to zero 
amplitude the initial data become, the solution can contain unbounded oscillations at some 
finite distance from the initial line; hence, the Cauchy problem is poorly posed, in gen- 
eral, for elliptic equations. This feature has been evidenced in numerical solutions 
of the blunt-body inverse problem (ref. 22, p. 452) as well as in the present problem. 


22 



When the method of lines is applied to the solution of Laplace's equation in a rec- 
tangular region one finds that the solution is composed of a complementary part which is 
a sum of exponential eigenfunctions, and a particular integral. The largest eigenvalue of 
the complementary part is proportional to the number of lines so that round-off errors 
can grow exponentially like 10 "i exp(Nx), where j is the number of decimal figures 
carried in the computation, N is the number of lines and x is the distance along lines 
from the initial data line, the lines being parallel to the X-axis. Hence, there exists an 
optimum number of lines which represents a compromise between decreasing discretiza- 
tion errors and increasing round-off error growth. Such a compromise is contained in 
most numerical procedures for the solution of differential equations. For example, in 
a standard solution by the grid relaxation method for Laplace's equation, the discretiza- 
tion error is 0(h2) whereas the round-off error is 0(h _ 2) (ref. 49, p. 482). Hence, the 
grid method has an algebraic round-off error growth whereas the method of lines has 
exponential error growth. 

The discussion is not intended to discourage the use of the method of lines; rather, 
it is intended to provide an understanding of difficulties if they arise in applications to 
nonlinear problems. One should expect to encounter instability where (1) the shock wave 
is a large distance from the body, a condition which occurs at low Mach number 
(Moo « 1.5) and on the leeside of the body at large incidence; (2) when cross-flow grad- 
ients are large and require a fine spacing of lines, which increases the eigenvalues; and 
(3) when a small number of figures is carried in the computation. These three condi- 
tions correspond to x large, N large, and j small, respectively, in the linear 
problem. 

In the conical flow problems treated so far, it appears that N = 19 is about the 
maximum number of lines which can be used without encountering sudden error growth. 
This observation applies to cases involving moderate -to -large Mach numbers (M^ ~ 2), 
angles of attack up to relative incidence of one, and a 60 -bit computer word length (about 
15 decimal figure accuracy). Usually, a much smaller number of lines (N ~ 9) provides 
sufficient accuracy for, say, circular cones, elliptic cones of moderate axis ratio 
(1.0 = a/b = 1.5), and delta wings at large Mach numbers (M*, 3). For investigations of 

cases which are more severe, for example, low Mach numbers, one will have to be con- 
tent with the accuracy provided by a smaller number of lines in order to avoid the large 
eigenvalues and thus suppress the error growth. 

Force Coefficients 

The force coefficients are computed by numerical integration of the surface pres- 
sures. The equations for these coefficients are given in appendix C. The reference area 
for the force coefficients is taken as the base area for the elliptic and circular cones and 
the plan area for wings. No force coefficients are presented in the results. 
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RESULTS AND DISCUSSION 


Circular Cone 

The circular cone is perhaps the most basic conical configuration and, as a con- 
sequence, has received more extensive treatment in the literature than other conical 
bodies. At small incidence it is one of the simplest conical flows to calculate and good 
experimental data are available to assess the validity of the computations. Moreover, at 
large values of incidence it exhibits the features which often complicate the computa- 
tion of many conical flows. An extensive tabulation of the flow about circular cones, not 
only surface quantities but the flow variables throughout the field as well, is presented 
in references 13 and 50. The latter results have been computed by the method of lines 
as presented in reference 2. Computations by the present method give comparable 
results. 

Tracy (ref. 51) presents results of an experimental study of a circular cone with 
semiapex angle of 10° at a Mach number of 7.95. Included are the measurements which 
delineate the shock structure and the viscous boundary; the region between the body and 
the viscous boundary is dominated by viscous effects. Figure 4 shows these measured 
boundaries and the shock shape computed by the method of lines for N = 15 at angles 



Figure 4 .- Cross section of shock for circular cone at incidence. = 7*95 J 
0 = 10°; N = 15. Experiment from reference 51. 
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Figure 4.- Concluded. 

of attack of 4° and 8°. It is seen that the computed shock shape is in good agreement 
with the measured results over the windward side but is less satisfactory on the leeward 
side where the viscous effects distort the flow. The experimental results show that the 
viscous region becomes more prominent as the angle of incidence increases and the com- 
puted shock wave on the lee side departs more from the measured values at the higher 
incidence. The corresponding surface pressure distributions are shown in figure 5 as a 
function of the circumferential angle <f> where <p is -90° in the windward plane of 
symmetry. Agreement with experiment is less satisfactory in regions where viscous 
effects are significant. Thus, the measured and calculated values show poorer agree- 
ment on the lee side and depart more at the larger angle of attack. 

Difficulties at large rela tive incidence .- As the angle of incidence is increased, 
conditions are encountered where it becomes difficult to obtain a converged solution. 

These conditions generally occur where a/Q 0 is near or exceeds unity. Three factors, 
singly or together, contribute to the difficulty in obtaining a solution: (1) the shock wave 
on the lee side approaches tangency with the stream Mach cone, (2) the cross flow becomes 
locally supersonic, (3) the computed pressure in the flow field or on the body becomes 
very small. It may be seen in figure 5 that the surface pressure on the lee side is 
tending toward the free -stream value Cp = 0 at a = 8°. Similarly, the shock -wave 
pressure coefficient (not shown) at the leeward symmetry plane is also near zero; hence, 
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Figure 5.- Circumferential pressure distributions on circular cone at incidence. 

Ho = 7-95; 9 = 10°; N = 15. 


that portion of the shock is approaching tangency with the free -stream Mach cone, which 
is, in turn, a characteristic surface. The result is that the expressions for the 
^-derivatives approach the indeterminate form (0/0) and the problem becomes very ill- 
conditioned. The particular sequence of runs for these calculations with 15 lines was 
made in increments of 1° in a for values of a up to 9° by using the previously 
described extrapolation routine for generating the initial values of ? 7 g for the pivotal 
computation. The routine made an attempt at a = 10° and encountered difficulty; as 
explained earlier, the 1° a -increment was halved, and successful convergence was 
achieved for a = 9.5°. Then it was necessary to halve the increment twice more to 
avoid difficulties, so that a = 9.625° was the next converged solution, and so on. Thus, 
the calculations became very sensitive for a ^ 9.5°; that is, a high degree of accuracy 
in the shock shape is required to satisfy the convergence criterion. The cross-flow 
Mach number M c at the shock surface is shown in figure 6 for a = 4°, 8°, and 9.5°. 

It may be seen that M c becomes 1 at the shock for a = 9.5° to add to the difficulty 
in obtaining a converged solution in this instance. Once difficulty is encountered, the 
computations can often be continued by reducing the number of lines. Thus, as previously 
discussed, there is generally some trade-off between the resolution and the range of the 
computations. The reduction of the number of lines from 15 to 13 permitted calculations 
for increments in a of 1° up to 10°, and by half -degree increments up to 11°. There- 
after, smaller increments were required. The cross -flow Mach numbers M c at both 
the body and shock are shown in figure 7 for the circular cone with 13 lines at a = 11°. 

It may be seen that a fairly extensive supersonic cross flow has developed, the largest 
values of the Mach number being at the body; at the shock the cross-flow Mach number 
does not substantially exceed one. 
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A pressure minimum on the circular cones generally occurs away from the plane 
of symmetry, and the pressure coefficient becomes negative at the larger values of 
incidence. The sensitivity of the calculations to the accuracy of the shock location 
coupled with the small values of the pressures often lead to failure of the method through 
computation of a negative value of the pressure under these circumstances. 

Extrapolation of surface pressures to large angles of incidence.- Even though the 
computations break down at some value of a near or beyond d Q , useful results on the 
windward side can be obtained for larger values of a. Extrapolation of the results has 
proven practical with the present program as also in the similar computations of Jones 
(ref. 2). In figure 8 extrapolated values of the surface pressure coefficients on the wind- 
ward side of the circular cone are compared with the experimental values measured by 
Tracy. The extrapolated values at a = 12°, 16°, 20°, and 24° computed with the use of 
the Lagrange formula from calculations by the method of lines at a = 6°, 8°, and 9.5° 
are in excellent agreement with experiment. 



Figure 8.- Comparison with experiment of large -incidence windward-side pressures 
extrapolated from calculations at a = 6 °, 8°, and 9-5°* ~ 7*95; 6 = 10°; 

N = 15. 
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Artificial hump on leeward side. - An alternate method of extending the calculations 
to large angles of attack is to alter the geometry on the lee side of the body to avoid some 
of the computational difficulties already discussed. The success of such an approach 
hinges on the fact that the flow on the windward side of the body is, to a large extent, 
independent of the flow on the lee side. Thus, the flow on the windward side of the spec- 
ified body can be computed by selecting the shape on the lee side in such a way that one 
source of the computational difficulties is avoided; that is, the pressures are raised on 
the lee side. It may be noted that a supersonic cross flow may still exist along the sides 
but does not appear to limit the computations. For a conical body that was semicircular 
on the windward side and semielliptic on the leeward side, labeled bielliptic herein, com- 
putations have been made for angles of incidence much larger than is possible for a cir- 
cular cone. The ratio of the major axis to the minor axis of the ellipse was taken suffi- 
ciently large to avoid computational problems and was sometimes taken as a linear 
function of the angle of attack for simplicity. Results of these computations at angles 
of attack of 12° and 16° for a stream Mach number of 7.95 and a semicone angle of 10° 
are shown in figure 9. Figure 9 displays the circular body, the experimentally deter- 
mined viscous boundary, and the shock contour (ref. 51). The body shape on the lee side 

Experimental shock 



Figure 9»- Comparison of shock shape computation for bielliptic cone with experi me nt 
(ref. 51) for circular cone at large incidences. = 7-95; 0 = 10 ° ; N = 15 . 
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Experimental shock 



Figure 9 -- Concluded. 

used for the computations is also shown as well as symbols for the computed shock shape. 
For both computations the calculated shock shape on the windward side is in excellent 
agreement with the experimental values. The reasonably good agreement between exper- 
imental and computed values of the shock location on the lee side for a = 12° is 
apparently due to a fortuituous choice of the axis ratio a/b of the ellipse (1.52) and was 
not selected essentially to duplicate viscous boundary. The axis ratio for the case at 
a = 16° was 2.52. The corresponding values of the surface pressures are shown in 
figure 10. The computed values are in good agreement with experiment on the windward 
side at a = 12° but show some deterioration at the more extreme case with a = 16°. 

It should be noted that the elliptic hump has held the leeward pressures at a "safe" 
positive level. 
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Figure 10.- Comparison of pressure computation for bielliptic cone 
with experiment for circular cone at large incidence. 


Entropy Layer and Vortical Singularities 

In the inviscid flow about a conical body, there exists a thin layer adjacent to the 
surface where the entropy gradients are large (in fact, unbounded). This so-called 
entropy or vortical layer was first recognized by Ferri (ref. 8) and has been studied 
analytically in many subsequent papers (refs. 9 to 11 and 52 as well as additional papers 
referenced therein). The pattern of the isentropes, or cross-flow streamlines, is shown 
in figure 11 as projected onto a plane Z = Constant, and illustrates the nature of the 
entropy layer for a circular cone at small incidence. In this case the streamline in the 
windward symmetry plane wets the entire cone surface, and it carries the maximum 
entropy. All other cross-flow streamlines rapidly approach the surface, form a thin 
layer of large (unbounded) entropy gradients, and converge in the leeward symmetry 
plane at a nodal point of the isentropes, where the entropy is multivalued. The nodal 
point is called the "vortical singularity." 
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o Nodal point 
• Saddle point 



Figure 11.- Cross-flow streamline pattern for circular cone at small 
relative incidence. M m = 5°5 0 O = a = 2°; N = 17 . 

Entropy layer. - As explained in an earlier section entitled "Extrapolation to 
Surface," it is not necessary to account directly for the entropy layer. Nevertheless, 
because the independent variable £ is continuous in the present method, it is possible 
to refine the integration near the surface and resolve some of the detail of this thin layer. 
Similarly, the "BVLR" method^ can be used to resolve such details close to the surface, 
since the implicit nature of the method allows a nearly unlimited refinement of the grid 
size normal to the surface; such computations have been presented in reference 53 for an 
elliptic cone at zero incidence. Using a semidiscrete method somewhat similar to the 
present one, Ndefo (ref. 3) has also obtained some entropy layer calculations for a circu- 
lar cone. 

The computation of the flow in the entropy layer with the method of lines is tedious 
because e must be very small in order to carry the integration very close to the sur- 
face. This requirement, in turn, demands that the shock shape be precisely determined. 
The necessary precision of the computation requires many iteration cycles and the cal- 
culation is generally feasible only with the use of the modified Newton iteration technique 
previously described. For the final computations presented here, e was 10-1°. 

4 After Babenko, Voskresenskii, Lyubimov, and Rusanov, the co-authors of 
reference 13. 
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Calculations within the entropy layer were made for a circular cone with a semi- 
cone angle of 10°, a = 2°, and M w = 5. The density distribution between the body and 
shock for 0 = 0 (along the side of the cone) is shown in figure 12. The density 



Figure 12.- Density distribution across shock layer and entropy layer 
of circular cone, = 5*0; 0 = 10°; a = 2°; $ = 0°. 

monotonically increases from the shock to a region near the body where, in the entropy 
layer, the density decreases. The inset in figure 12 expands the abscissa by a factor 
of 50 and the ordinate by a factor of 10 to display the variation near the conical surface. 

The results for several values of 0 are presented in figure 13 where the quan- 
tity p - p 0 (p 0 is the surface value of the density) is presented as a function of £ on 
a log-log plot. Those results show that deep in the entropy layer, the density variation 
is of the form 

pw>,o = p 0 (<m + com" 

From figure 13, v ~ 0.19. Dr. R. E. Melnik has supplied the authors with an estimate 
of v based on his asymptotic analysis of the entropy layer (ref. 9). His estimate for 
this case is v = 0.22 ± 0.04 which is consistent with the present numerical results. The 
possible error indicated is a result of neglecting terms which are i n the 

estimate. 
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Figure 13-- Density variation in the entropy layer, = 5*0; 9 = 10°; a = 2°. 

Vortical singularity lift-off. - Ferri (ref. 8) suggested the possibility that the nodal 
singularity could lift off the conical body at large angles of incidence. Several authors 
have attempted to develop analytic solutions in the neighborhood of the nodal point; 

Melnik (ref. 52) and others have concluded from their analyses that the vortical singu- 
larity can, in fact, lift off the surface. 

Some effort has been expended to verify the lift-off phenomenon by exact numerical 
solutions, with questionable results. Gonidou (ref. 16) has attempted such calculations 
with the BVLR method, but his results indicated numerical instability, and hence were in 
doubt. It has been shown (ref. 53) that the conditions for lift-off are just those which 
cause the BVLR method to be unstable. Jones (ref. 2) has made computations which 
indicate that the vortical singularity was off the conical surface. His calculations are 
for a circular cone with a semivertex angle of 12.5° at a stream Mach number of 1.797. 
Jones presents computations for a/0 o up to 1.4 and his results indicate that lift-off 
occurs for a/Q 0 between 1.1 and 1.2. His integration step size was one-tenth of the 
shock -layer thickness and extrapolation to the surface was made from the last step. 

Thus, he extrapolated across the nodal point for the angle -of -incidence conditions where 
the vortical singularity is off the surface. 
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Computations have been made with the method of lines presented herein in an 
attempt to verify Jones' results. The calculations were made for a/B 0 up to 1.2. 

Sonic velocity was attained at the shock in the lee plane of symmetry for this condition 
and the computations could not be carried to larger values of a because of the sensi- 
tivity of the calculations to the shock shape. One set of computations was made with an 
integration step size of 0.1 with extrapolation to the surface from 0.1 to nearly duplicate 
the calculations of Jones. A second set carried the integration to a value of 0.05 with 
extrapolation from that value of £ to the surface. Both calculations indicated that lift- 
off had just barely occurred for a/0 o = 1.2. The results of the latter computations are 
shown in figure 14 where dv/d£ and v in the leeward plane of symmetry are plotted 
against £. The values of dv/d? were extrapolated (shown dashed) to £ = 0 with the 
use of the Lagrange formula and the corresponding values of v were found by integra- 
tion of dv/d£. The incipient lift-off condition occurs for = 0 at £ = 0; at angles 

of incidence beyond the lift-off condition, > 0 at the surface and a zero value of v 

occurs off the body (nodal point) as well as on it. Precise lift-off computations are 



z z 

Figure Ik.- Extrapolation of the normal velocity component and its derivative in the 
lee plane of sy mm etry of a circular cone near the condition for vortical singularity 

lift-off. = 1-797; e = 12.5 0 ; = 15°; n = 13. 
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difficult to make since the shock must be located with great accuracy in the region near 
the leeward plane of symmetry. In practice, computations are made for several values 
of ot near or beyond the incipient lift-off condition. In cases where the nodal singu- 
larity is off the surface, the quantity dv/ d£ is extrapolated across the singular point. 

The calculation of the flow beyond the incipient lift-off condition in this manner is 
questionable. The extrapolation of the derivative dv/d£ across the singular point in 
the leeward symmetry plane is one facet open to criticism since the discontinuity in the 
entropy across the singularity is not taken into consideration. That is, the entropy should 
be that of the leeward symmetry plane between the shock and the nodal point and that of 
the windward symmetry plane between the body and the nodal point. Likewise, there is a 
discontinuity in both p and u across the nodal point. From the second of equations (2), 
in the plane of symmetry (w = 0) and in the neighborhood of the nodal point (v — 0), the 
expression for the derivative becomes 


I ■ -*V 


Thus, dv/d£ is discontinuous like u at the nodal point and extrapolation of dv/d£ 
across it is invalid. 

The method of lines could be adapted to account properly for the discontinuity in the 
variables across the vortical singularity although it is not known whether it would indeed 
be successful because of the computational difficulties that are encountered at the large 
incidences where lift-off is supposed to occur. The computations would require extrapo- 
lation from some integration step to a value of £ at which the v component of velocity 
vanishes in the leeward symmetry plane, that is, to the location of the nodal point. The 
jump in entropy across the singularity is known and the jump in p and u can be com- 
puted. The integration could then be continued with the new initial values and finally 
extrapolated to the surface in the normal fashion. No attempt has been made to carry 
out such a sequence of computations. In view of the incorrectness of the extrapolation 
across the nodal point which has been used in the method of lines, and the instability of 
the BVLR method when incipient lift-off conditions occur, it is felt that the existence of 
a nodal point off the surface has not yet been conclusively demonstrated by numerical 
computation. 


Elliptic Cone 

Calculations have been made by the method of lines for the supersonic flow over an 
elliptic cone to compare with the experimental results of Chapkis (ref. 54). The cone 
angle d Q in the vertical plane of symmetry is 5.97°, the axis ratio a/b is 2, and the 
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stream Mach number is 5.8. The computed and experimental surface pressure coeffi- 
cients are shown in figure 15 for a of 0°, 2°, 4°, and 6°. The agreement between the 



Figure 15.- Computed, and experimental circumferential pressure distributions for 
elliptic cone at incidence. M m = ^. 8 ; 0 Q =» 6°; ^ = 2; N = 19* 

computations and experiment is generally good but is less satisfactory on the leeward 
side where the viscous buildup raises the level of the measured pressures. Both the 
computed and experimental results show a relatively constant value of the pressure on 
the windward side and a rapid expansion around the side of the elliptic body, a supersonic 
cross flow occurring for the case at an angle of attack of 6°. Typical of the inviscid cal- 
culations for both the circular and elliptic cones is the pressure minimum which occurs 
away from the horizontal plane of symmetry when the relative incidence a/ 9 Q 
approaches unity. Convergence of the solution for angles of incidence beyond 6° 

becomes difficult, similar to the case of a circular cone for ^- = 1, since the imbedded 

u o 

regions of supersonic cross flow and low pressures are encountered. 

Extrapolated surface pressures. - Extrapolated values of the surface pressure coef- 
ficient on the windward side of the elliptic cone, using the Cp values for a = 2°, 4°, 
and 6°, are compared with the experimental results of Chapkis in figure 16. The agree- 
ment is excellent at a = 8° and 10° but shows some deterioration at a = 14°. 
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Figure 1 6.- Comparison with experiment of large incidence windward-side pressures 
extrapolated from calculations at a = 2°, 4°, and 6°. = 5-8; 9 q « 6°; 

- = 2; N = 19- 
b 

Cross-flow streamline pattern. - In figure 17 the cross-flow streamlines together 
with the cross section of the body and shock are shown for the elliptic cone at angles of 
incidence of 2° and 6°. At an angle of attack of 6° the relative incidence a/6 0 is 
about 1. The streamline contours, which correspond to constant values of the entropy, 
were obtained by linear interpolation. Those contours either near the windward saddle 
point or the leeward nodal point are difficult to obtain with precision and are shown 
dashed in some cases where their location is uncertain. It may be seen that a nodal 
point (shown with an open symbol) occurs in the windward and leeward planes of sym- 
metry and a saddle point singularity (shown with a filled symbol) on the windward side. 
At zero incidence the saddle point lies on the major axis and moves to the windward side 
with increasing incidence. The cross -flow streamline patterns are similar on the wind- 
ward side for the two cases; however, the general character of the stream lines on the 
leeward side changes as the relative incidence is increased. 
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O Nodal point 
• Saddle point 




Figure 17.- Cross-flow streamline pattern for elliptic cone at incidence. 

H. = 5 . 8 ; e 0 ~ 6 °; | = 2; N = 19- 


Computation history.- The elliptic cone calculations at a = 0° were made by 
starting from a circular cone (a/b = 1) with semiangle equal to the desired semiangle 
(11.8°) in the plane of the major axis of the ellipse. The axis ratio a/b was changed 
in increments until the desired configuration was obtained. The solutions for various 
angles of attack were then computed, beginning with zero angle of attack. Table I gives 
the history of the automatic increases in a together with the number of Newton iteration 
cycles (N + 1 integrations) and the number of integrations of the coupled 5N differential 
equations required to converge the shock shape for each a by using the regular Newton 
iteration. The convergence criterion for the normal velocity component was set at 
e = 10 and the final value of max|v-(0)j is also shown. Each integration from the 

shock to the body was made in 12 integration steps with a step size of 0.1 to a value of £ 
of 0.1, and steps of 0.05 and 0.025 thereafter. The total central processing time, 
excluding compilation, was 740 seconds (Control Data 6600 computer). Table it gives the 
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TABLE I.- COMPUTATION HISTORY FOR ELLIPTIC CONE 


WITH REGULAR NEWTON ITERATION 
[n=17; M 0O = 5.8; ^=2; d Q = 6°] 


Angle of 
incidence* 
a, deg 

Number of 
Newton cycles 

Number of 
integrations 

maxh 
i ' 

/i(0)| 

0 

0 

i 

7.7 

X 

10-4 

.05 

1 

19 

3.7 

X 

10-5 

1 

1 

19 

5.7 

X 

10-4 

2 

2 

37 

4.9 

X 

h- * 

o 

1 

CJ1 

3 

1 

19 

3.8 

X 

10-5 

4 

1 

19 

2.6 

X 

10-5 

5 

1 

19 

1.6 

X 

O 

1 

6 

2 

37 

1.7 

X 

10-5 

Total . . . 


170 





TABLE H.- COMPUTATION HISTORY FOR ELLIPTIC CONE 
WITH MODIFIED NEWTON PROCEDURE 
[N=17; M 0O = 5.8; g = 2; & Q = 6°] 


Angle of 
incidence* 
a, deg 

Number of 
Newton cycles 

Number of 
modified 
Newton cycles 

Number of 
integrations 

max|v^(0)| 

0 

0 

0 

1 

7.7 X 10-4 

.05 

1 

0 

19 

3.7 X 10-5 

1 

1 

0 

19 

5.7 X 10-4 

2 

1 

2 

21 

2.6 X 10-4 

3 

1 

0 

19 

4.2 X 10-5 

4 

1 

0 

19 

4.2 X 10-5 

5 

1 

0 

19 

1.6 X 10-4 

6 

1 

1 

i 

20 

2.9 X 10-4 

Total . . . 



137 
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computation history for the same set of computations except that the modified Newton 
procedure was employed; that is, the partial derivatives employed in the Newton iteration 
were not recomputed in cases where more than one iteration was required. The central 
processing time in this case was 595 seconds. Comparison of tables I and II shows that 
worthwhile savings in computing can be obtained with the use of the modified Newton 
method. The superiority of this method is more dramatic for sequences of computations 
where the regular Newton method requires two or more cycles for each value of the 
parameter. Tables I and II show an increase in the number of integrations for angles of 
attack of 2° and 6°. The increase at an angle of attack of 2° is due to the fact that the 
initial values of ri s are found by linear extrapolation of those for a = 0° and a = 1° 
whereas the initial values of t? s at larger angles of incidence are obtained by quadratic 
extrapolation and consequently provide a better estimate for the pivotal computation. At 
the angle of attack of 6°, ct/0 Q is approximately one, a supersonic cross flow has 
developed and convergence, as already discussed for the case of the circular cone, is 
more difficult to achieve. 

Effect of number of lines. - In problems involving severe cross-flow gradients, it 
would be desirable to use a relatively large number of lines to resolve the detail of a 
solution. However, the maximum number that can be used profitably in a given compu- 
tation is restricted by the computing time and by instabilities that enter the solution, as 
explained in an earlier section. A maximum of about 19 lines has proven feasible for the 
method of lines in the form presented herein. 

Figures 18 and 19 illustrate the changes in the solution with a change in the number 
of lines for the elliptic cone with a/b = 2, a = 2°, = 5.97°, and a stream Mach num- 

ber of 5.8. The surface pressure distribution is shown in figure 18 for N = 5, 9, and 17 
and the corresponding shock shapes, together with the body shape, are shown in figure 19. 
It may be seen that the results for N = 5 and N = 9 oscillate about those for N = 17. 
The computations for N = 5 clearly establish the trend of the solution whereas those 
for N = 9 are in very good agreement with the 17 -line case. The central processing 
time, excluding compilation, for the same sequence of computations beginning with a 
circular cone was 46, 103, and 322 seconds for N = 5, 9, and 17, respectively. 

Elliptic cone with large axis ratio. - The pressure gradients in the cross-flow direc- 
tion for the elliptic cone become large as a/b becomes large; thus, a relatively large 
number of lines to construct the solution to good accuracy are dictated. However, the 
instabilities which arise with a large number of lines force the use of fewer lines than 
would otherwise be desirable. Computations for an elliptic cone with N = 17 at zero 
incidence are compared with the experimental results of Martellucci (ref. 55) in figs. 20 
and 21; the semicone angle in the plane of the major axis is 30.05°, a/b = 3.68, and 

= 3.09. The general agreement of the computed and measured pressures (fig. 20) is 
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good but the computed results clearly exhibit some roughness. The roughness of the 
shock shape (fig. 21), and, in particular, the too large value of the shock coordinate in the 
plane of symmetry of the major axis is typical of cases where too few lines have been 
employed. (See fig. 19.) It is possible to continue the solutions to larger values of a/b 
but converged solutions become more difficult to obtain, the "roughness" of the solution 
becomes more pronounced, and program failure ultimately follows. 

Comparison of an elliptic cone computati on with other methods.- The saddle and 
nodal point singularities, as seen in figure 17, arise when the w-component of velocity on 
the body vanishes. This component of velocity, which is one of the three velocity com- 
ponents of the coordinates used herein, is normal to a ray from the apex. Babenko 
(ref. 53) employed a cylindrical coordinate system to compute the flow about an elliptic 
cone and erroneously concluded that a singularity occurred where the circular -cylindrical 
component of velocity w c vanished and thereby computed a discontinuity in the other 
two cylindrical velocity components. Jones (ref. 2) similarly used a circular -cylindrical 
coordinate system but calculated the solution correctly. The semicone angle in the plane 
of the major axis is 35°, a/b =2, a = 0, and = 10 for this case. The circular - 
cylindrical velocity components computed by Babenko, Jones, and by the present method 
are shown in figure 22. The results of Jones and the present calculations are in 
agreement. 

Elliptic cone at yaw. - The designation "elliptic cone at yaw" refers to one for 
which the stream velocity vector is in a plane parallel to the major axis; p is the yaw 
angle for this configuration. Chapkis (ref. 54) has made measurements of the surface 
pressures for the yawed cone. The conical body is the same one for which the angle-of- 
attack measurements were made; the semicone angle in the plane of the minor axis 
is 5.97°, a/b = 2, and the stream Mach number is 5.8. The computed surface pressure 
coefficients are compared with the experimental results in figure 23 for angles of yaw 
of 2°, 4°, 6°, and 8°. The fairing of the experimental data was taken from reference 54. 
Note that this faired data does not have a zero slope at <p = 0° and 180° as should be 
the case. The pattern of the flow is somewhat different from that of the elliptic cone at 
incidence. The w component of velocity for the 4° yaw condition is shown in figure 24. 
The nodal point, which is located where the w component of velocity vanishes outside 
the plane of symmetry, is located on the leeward side of the conical body. The cross - 
flow streamlines together with the body and shock shape are shown in figure 25 for yaw 
angles of 2°, 4°, 6°, and 8°. The streamline contours near the nodal singularity are 
difficult to establish accurately and consequently are shown dashed. The saddle points 
are shown with a filled symbol and the nodal point with an open symbol in the figures. At 
zero yaw the nodal point is located in the plane of the minor axis and as the yaw angle is 
increased, it moves toward the leeward plane of symmetry. The surface entropy is dis- 
continuous across the nodal singularity; the surface entropy in the region from the 
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windward saddle point to the nodal point (w > 0) corresponds to the value of the entropy 
in the windward plane of symmetry whereas that from the nodal point to the leeward 
plane of symmetry (w < 0) corresponds to the value in the leeward plane of symmetry. 
Along with the discontinuity in surface entropy is a discontinuity in the u-component of 
velocity and the density. The surface distribution of the u-component of velocity for the 
4° yaw is shown in figure 26. As the angle of yaw increases, the nodal point continues 
to move toward the leeward side and at some yaw condition apparently moves to the lee- 
ward plane of symmetry. At an angle of yaw of 8° (fig. 25), the nodal point appears to 
have coalesced with the leeward saddle point; however, it is impossible to obtain any 
resolution of the transition from the present computations. 



Figure 26 .- Circumferential distribution of u-component of velocity at surface 
of yawed elliptic cone, illustrating discontinuity at nodal point. = 5-8; 

9 0 = 11.8°; | = 2; \|c = 4 °; N = 15. 

Conical Delta Wings 

The present method has been applied to the problem of conical delta wings with the 
shock wave attached at the sharp leading edges. In this problem the compression surface 
is independent of the other surface, and the two can be treated separately. The method 
of lines, as formulated herein, is only applicable to the compression side where the 
shock wave forms the outer boundary. Nevertheless, the compression -side results are 
valuable from both theoretical and practical standpoints. As mentioned earlier, the 
results from the present method provide a check for more complex three-dimensional 
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calculation methods, and it will be seen that they corroborate results from other conical 
methods and experiment, too. Practically speaking, at large supersonic Mach numbers 
nearly all the aerodynamic -force contribution comes from the compression side. 

In these problems the cross-flow component Vv^ +w^ is supersonic at the 
leading edge and remains supersonic for some distance toward the wing center line. The 
partial differential equations governing the conical flow are of the hyperbolic type in 
regions where the cross flow (that is, the velocity component on the spherical surface) 
is supersonic, and a two-dimensional conical method of characteristics can be used to 
construct the exact flow in such regions. Maslen (ref. 28) was the first to describe this 
approach for conical delta wings, and it has been exploited extensively in reference 41 
for the flow about elliptic cones at large angles of attack. A conical method of character- 
istics has been developed by Chiang and Wagner (ref. 56); they supplied the conical char- 
acteristic results presented in this paper. 

Vincenti and Fisher (ref. 57) proposed an approximate method which can also be 
used in regions of supersonic cross flow; it is analogous to the familiar shock-expansion 
method for two-dimensional and axisymmetric flows. The approximation reduces to a 
pair of nonlinear ordinary differential equations that are numerically integrated along the 
surface in the spanwise direction from the leading edge inward to the cross -flow sonic 
point. Results obtained by Richard D. Wagner, Jr., by that method are presented and 
are referred to as a "conical shock expansion" method. 

The supersonic cross flow in the delta -wing problems presents no difficulties for 
the method of lines; this result is contrary to the case of the circular or elliptic cone at 
incidence where convergence becomes more difficult when a region of supersonic cross 
flow occurs. 

Parabolic-arc cross section . - Reference 17 presents calculated results for a num- 
ber of conical delta wings obtained by the "BVLR" method. Computations for delta wings 
with both flat and parabolic -arc cross sections presented in reference 17 are compared 
herein with results from the present method and other methods. Planform and section 
views of a thin parabolic cross-section wing are shown in figure 27 with the calculated 
shock shape for a = 10°. The conditions for this problem are M m = 4, A = 50°, 

9 q = 3°, and N = 12. The cross-flow sonic line is shown as a heavy dashed line in the 
section view. 

Equal line spacing was used for the wing computations with the method of lines; 
that is, equal increments were used around the wing contour on the surface r = 1. 

Thus, when projected onto the plane Z = 1 in which the body cross section is prescribed, 
the lines appear to spread as the leading edge is approached. The solid curve for the 
shock (fig. 27) represents Voskresenskii's results, and the squares are the present 
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Voskresenskii (ref. 17) 

□ Method of lines, IM * 12 
O Conical method of characteristics (ref. 56) 



Figure 2J . - Comparison of computations for shock shape and spanwise pressure distribution 
for a conical delta wing. Parabolic-arc cross section; = 4°; 0 = 3°; A = 50°; 

N = 12. 

results for N = 12. The replotting of Voskresenskii' s shock shape is probably not as 
accurate as that for the pressures because although Voskresenskii' s shock appears to lie 
slightly inside the present result near the center line, shock pressures by both methods 
(not shown here) agree very well except near the leading edge. The conical method of 
characteristics also gives a shock shape from the leading edge to the cross-flow sonic 
lines, and the results (not shown) are in agreement with the present results. 

The surface pressure coefficient is plotted against x, a nondimensional spanwise 
coordinate; x = 0 at the wing center line and x = 1 at the leading edge. The vertical 
dashed lines indicate the spanwise location of the cross-flow sonic point^ at the surface 
and, hence, the limit of applicability of both the conical method of characteristics (MOC) 
and the conical shock -expansion method. The method of lines, conical characteristics, 
and conical shock expansion agree very well in the supersonic cross -flow region. 
Voskresenskii' s results are generally quite close to the present results. However, his 
results for a = 5° are generally higher than the other computations in the outboard 

5 As the angle of attack is increased, the cross -flow sonic point moves toward the 
leading edge and reaches it slightly before leading-edge detachment occurs. 
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region except at the leading edge. At the leading edge, where the pressure and other 
variables can be calculated accurately from the algebraic shock conditions, results from 
all theories should agree exactly yet Voskresenskii's result for the pressure is noticeably 
low. His pressure calculation just inboard of the leading edge rises slightly above that 
of the three conical methods, but then appears to agree well with the present results in 
the central subsonic cross-flow region. It is not entirely clear whether the generally 
small discrepancies in the supersonic cross -flow region are due to the error in the 
leading-edge boundary condition. Of particular note is the very good agreement of the 
comparatively simple conical shock -expansion calculations with the other results. 


Circular -arc cross section. - Reference 58 presents experimental data for various 
conical delta wings with attached leading-edge shocks; one of the models had a circular- 
arc cross section, a good test problem for the various computational methods. Figure 28 
shows the shock shape computed by the method of lines for N = 12 and spanwise pres- 
sure distributions compared with experimental and other computed results. The wing 
thickness and the Mach number are both about twice that of the computations in figure 27. 
It can be seen that the cross-flow sonic line is closer to the wing center line; this result 
is due largely to the higher Mach number. 



3-dimensional method of characteristics (ref. 18) 

□ Method of lines, N = 12 
O Conical method of characteristics (ref. 36) 

• Conical shock expansion 

Experiment (ref. 58) 



Figure 28.- Comparison of computations and experiment for shock shape and spanwise 
pressure distribution for conical delta wing. Circular-arc cross section: 

M ro = 8.1; 9 = 6.54°; A = 50°; N = 12. 
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Hatched bands are used to represent the measurements of reference 58 for several 
different Reynolds numbers and stations downstream from the wing apex. The measured 
pressures are higher than the inviscid predictions everywhere on the surface; however, 
the agreement is generally good except near the leading edge. The higher pressures are 
caused by the hypersonic boundary -layer displacement effect which is most pronounced 
near the leading edge (x = 1). Computations by a three-dimensional method of character- 
istics described in reference 18, conical method of characteristics, and conical shock 
expansion are shown for a = 10°. The agreement is generally good in the supersonic 
cross -flow region; in the region of subsonic cross flow, three-dimensional method of 
characteristic calculations are generally lower than those computed by the method of 
lines. At the smaller angles of attack, the three-dimensional method of characteristics 
and method of lines computations are in excellent agreement and lie at the lower limit of 
the experimental data. 

The cross-flow streamline pattern for the circular-arc delta wing is shown in fig- 
ure 29. The streamlines show very little curvature; however, the computations were 



Figure 29*- Cross-flow streamline pattern for conical delta wing. 

Circular-arc cross section; M OT = 8.1; 0 = 6.5^°; A = 50°; 

a = 10°; N = 12. 

made by linear interpolation and consequently the finer detail may be lost, particularly in 
the region of the nodal point. The nodal point (shown with an open symbol) lies in the 
plane of symmetry and the surface entropy is equal to that at the wing leading edge. 

Spanwise pressure distributions for the circular-arc delta wing for = 5.08 are 
shown in figure 30. Both the three-dimensional method of characteristic calculations and 
those by the method of lines (N = 10) are in good agreement with the experimental results 
of reference 58 over most of the wing and the agreement near the leading edge is better 
than that shown in figure 28 since the viscous interaction effects are not as severe at the 
lower Mach number. The shock shapes for a =4.5° and a = -5° are also shown in 
figure 30. The three-dimensional method of characteristic computations are in good 
agreement with those by the method of lines with some minor differences near the plane 
of symmetry. 

Flat delta wing . - The classical problem of the flat -plate delta wing has several 
interesting features. There are three distinct regions to the flow field which exhibit dif- 
ferent flow characteristics. The region bounded by the wing surface, the cross -flow 
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3-dimensional method of characteristics (ref. 18) 

□ Method of lines; N - 10 
\SS^ Experiment (ref. 58) 



Figure 30.- Comparison of computations and experiment for shock shape and spanwise 
pressure distribution for conical delta wing. Circular-arc cross section; 

M m = 5.08; 9 q = 6.54°; A = 50°; N = 12. 

sonic line, and shock is one in which all the flow quantities p, p, V, and S are con- 
stant since the shock wave is planar from the leading edge to the cross -flow sonic line. 
The isentrope extending from the intersection of the cross-flow sonic line with the shock 
is a boundary of the second region. The region bounded by this dividing streamline, the 
cross -flow sonic line, and the wing surface is one of varying flow properties but is isen- 
tropic. The third region bounded by the dividing streamline, the shock, and the plane of 
symmetry is one of variable entropy. The exact inviscid pressure distribution should 
exhibit a "corner" of discontinuous slope at the cross-flow sonic point. These weak 
singularities are not treated in any special way in the present method and the same is 
assumed to be true of reference 17. 

The computations for the flat -plate conical delta wing shown in figure 31 are for 
A = 50° and = 4. The method of lines computations and those of Voskresenskii 
produce nearly identical results. The open circles are used to show the actual location 
of the computational lines in the present method for N = 12; the results are surprisingly 
smooth around the sonic line. Although the assumed starting shock shape (eq. (7)) is 
analytic, the final converged shock points lie very close to a straight line between the 
leading edge and the cross -flow sonic line as they should. A comparison is made in 
reference 20 of the pressure distribution and the shock shape for this wing computed by 
the method of Kutler and Lomax and with the method of lines; the agreement is 
excellent. 
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— Voskresenskii (ref. 17) 

O Method of li nes, N = 12 



Figure 31.- Comparison of computations for shock shape and spanwise pressure 
distribution for conical delta wing. Flat-plate cross section; M OT = 4; 
0 O =0°; A = 50°; N = 12. 


In the methods of references 18 and 30, the flow properties are set constant in the 
supersonic cross flow and in reference 30; an effort was made to account for the weak 
singularities. Even so, the latter results (solid circles) are in considerable disagree- 
ment with the other methods, particularly at the larger values of a. It is also note- 
worthy that although Babaev (ref. 30) attempted to account for the singularities, his 
results for the surface pressure distribution appear to be very smooth at the sonic point, 
without a corner. Results by the method of reference 18 (three-dimensional method of 
characteristics) are shown in the Cp plot for a = 15°, and they are in general agree- 
ment with the present results and those of reference 17; apparently, the results of refer- 
ence 30 are erroneous. 

The cross-flow streamline pattern for the flat -plate delta wing at an angle of attack 
of 15° is shown in figure 32. The cross-flow streamline emanating from the intersection 
of the sonic line and the shock is the dividing streamline which separates the outboard 
constant entropy region from that of variable entropy. No attempt was made to trace the 
cross -flow streamlines in the constant -entropy region. 
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Figure 32.- Cross-flow streainline pattern for conical delta wing. Flat-plate cross section; 

^=8.1; 0 O = 0°; a = 15°; A = 50°; N = 12. 


The method of lines computations for a flat conical delta wing are compared with 
experimental values of reference 58 for = 5.08, A = 50°, and a = 14° in figure 33. 
The experimental values (shown hatched) are measurements from several streamwise 
stations. The method of lines computations generally fall at the lower limit of the exper- 
imental data. 

The flow over a flat delta wing, = 4, A = 50°, has been computed for large 
angles of attack approaching the condition of shock detachment from the wing leading 
edge. Both the shock and surface pressures are shown in figure 34 for a = 15°, 20°, 
and 22.4°. The subsonic cross-flow region increases in extent as the angle of incidence 
increases and just before the shock detachment condition, the cross flow is everywhere 
subsonic. The shock detachment condition was attained for a slightly greater than 

M? Experiment (ref. 58) 

O Method of lines, N = 10 



Figure 53-- Comparison with experiment for spanwise pressure distribution for conical 
delta wing. Flat-plate cross section; M = 5.08; 0 = 0°; A = 50°; a = 14°; 

N = 10. 
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— O- Body 
— O— Shock 



Figure 34.- Calculated spanwise pressure distribution on flat delta wing for conditions 
approaching leading-edge shock detachment. M ro = 4°; 0 Q = 0°; A = 50° j N = 12. 

22.4°. It may be seen that even when the cross flow is all subsonic, there is little varia- 
tion of the pressure between the shock and wing surface at a given spanwise station. 

Wing with reverse curvature. - Measurements of surface pressures are presented 
in reference 58 for a delta wing with a bell-shaped cross section that has a change in sign 
of the curvature. The cross-section shape is given by 

y 0 = -tan 0 o (l - £2 2 ) 3 
x o 

£2 = — 2 tan 9 
b 

where b determines the spanwise location at which y Q and its first two derivatives 
vanish. The results presented here are for the configuration designated (B-W)j in refer- 
ence 58 with b = 1; thus, the slope goes to zero at the wing leading edge. The shock 
shape and pressure distributions are shown in figure 35. The conditions are M TO = 8.1, 

A = 50°, a = 4°, with 9 Q = 4.8° and 9 Q = 12.5°. The computations for 9 Q = 4.8° with 
N = 12 show a smooth shock shape and pressure distribution. The computations for 
& Q = 12.5° with N = 9 show some "roughness” in both the shock shape and the pressure 
distribution. The computed values of the shock coordinates for 9 Q = 12.5° are connected 
with straight lines to show more clearly the "roughness." The computed pressures follow 
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M? Experiment (ref. 58) 
O N = 9 



0 .2 .4 .6 .8 1.0 


x 

(a) Comparison of computations and experiment for spanwise pressure distribution. 



Figure 35-- Effect of increasing thickness on flow over a conical delta wing 
with reverse cross-section curvature. M m = 8.1; A = 50°; a = 4-°. 


the general distribution of the experimental data which is shown hatched. Noteworthy* is 
the very abrupt thickening of the shock layer in the region where large pressure gradients 
occur. Mead and Koch (ref. 58) suggested that a shock occurs in the cross flow. This 
possibility seems likely, in which case the method of lines is unsuitable for computing 
this flow. Attempts to compute the d Q = 12.5° case for N > 9 failed since the compu- 
tations showed increased roughness and convergence was not attained. The body- 
oriented coordinate system used here is not well adapted to configurations with reversed 
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curvature because the lines tend to converge as the shock is approached, with a resultant 
increase in computational difficulties. The question of an internal shock could probably 
be settled by tracing the characteristics in a cross-flow MOC computation. Kutler's 
method (refs. 19 and 20) would likewise be suitable for study of this case. 

Convergence history.- Figure 36 illustrates the convergence history for the flat- 
plate delta wing at a = 15° shown in figure 31. The surface pressure distribution 



Figure 36.- Surface pressure convergence history for flat delta wing. 

14^=4°; 9 0 = 0°; . A = 50°; a = 15°; N = 12. 

obtained with each new pivotal shock shape, including the initial approximation given by 
equation (7), is shown. The key of the figure gives the maximum normal velocity magni- 
tude, at the end of each cycle and a designation indicating whether the cycle was a regular 
or modified Newton cycle. 

Convergence with increas ing N .- Figure 37 provides an illustration of the accuracy 
that can be achieved with small values of N. A comparison of the surface pressures and 
the cylindrical cross -flow velocity component w c with N = 2, 4, and 8 is made with 
the computations of Voskresenskii for the flat delta wing at a = 15°. The pressure is a 
quantity that requires no correction after extrapolation to the surface whereas the veloc- 
ity w c is a surface -corrected quantity. With only one line placed between the leading 
edge and the wing center line (N = 2), the accuracy is surprisingly good. More lines are 
used, however, to give better resolution. 
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X 

(a) Comparison of computations for spanwise pressure distribution. 



o .2 .4 .6 .8 1.0 


X 

(b) Comparison of computations for spanwise distribution of cylindrical 

velocity component w c < 

Figure 37.- Convergence with increasing N for flat delta wing. M = 4 



Variable Line Spacing 


The computations presented herein were made for an equal line spacing in the arc 
length coordinate r; that is, the computational coordinate £ was equal to r. It would 
appear desirable to stretch the computational coordinate £ compared with the arc 
length t to bunch the lines in the region where the gradients are large and the coordi- 
nate transformation £ = |(r) was included in the program for this purpose. When the 
geometry and computational lines are projected onto a plane Z = Constant to appear 
like a cylindrical coordinate system, the line spacing is nonuniform in the polar angle <p. 
The smallest values of A 0 on an elliptic cone do, in fact, occur in the region of maxi- 
mum surface curvature and flow gradients, at the extremes of the major axis. The pre- 
viously described instabilities that arise with very small line spacing impose a limit on 
the amount of local "bunching" of the lines which can be utilized. Since about 13 to 17 
lines are required, in general, for elliptic cones, further bunching of the lines in the large 
curvature region can push the calculation beyond the stability limitation on At. Within 
this limitation some calculations have been made for elliptic cones in which the £ coor- 
dinate was stretched in the region of the major axis; however, the results were inconclu- 
sive regarding the merit of using such a nonuniform line spacing in the present "natural" 
coordinates. 


CONCLUDING REMARKS 

The method of lines has proved to be a useful and versatile procedure for con- 
structing the numerical solution to conical flow problems. It has the merit of a relatively 
simple mathematical formulation and can be adapted to compute the flow about a variety 
of conical bodies. The method has been applied to circular and elliptic cones at incidence 
and to the compression side of several conical delta wings with an attached leading-edge 
shock. Numerous comparisons have been made with experiment which serve to validate 
the method. 

Several other approaches for calculating conical flows have been compared with the 
method of lines and were found to be in agreement. Two of these approaches are three- 
dimensional computations that attain the conical solution asymptotically, whereas two 
others are strictly two-dimensional conical methods. The conical flow problems provide 
excellent check cases for the more general three-dimensional computational methods 
because a conical flow is three dimensional but can be solved by two-dimensional methods. 

Computations have been made for a circular cone to compare with experimentally 
determined values of both the shock shape and the surface pressures. The computed 
values are in good agreement with the measured results over the windward side but 
agreement is less satisfactory on the leeward side where a viscous buildup occurs to 
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raise the level of the measured pressures. The region dominated by the viscous buildup 
increases in extent with increasing incidence and the computed results for the leeward 
side show poorer agreement with experiment at the larger angles of attack. Comparison 
of computed and measured pressures for an elliptic cone of axis ratio a/b = 2 similarly 
show good agreement over the windward side. 

Whenever a supersonic cross flow occurs around circular or elliptic cones, con- 
vergence of the normal components of velocity at the body surface to a zero value is 
more difficult at lines in that region. Convergence is possible, however, if the super- 
sonic cross -flow region is not extensive. Low pressures, supersonic cross flow, and 
large cross -flow gradients are factors which limit the range of parameters for which 
computations can be made. These factors tend to limit the angle of attack to values such 
that the leeward side is not appreciably in wind shadow, the stream Mach number gener- 
ally greater than 2 but dependent upon cone angle, and the axis ratio a/b for the elliptic 
cone generally less than approximately 3. These limits are very approximate and depend 
too on the convergence criteria one is willing to accept. 

Cross -flow streamline patterns for the elliptic cone both at incidence and yaw are 
presented to illustrate the change in the flow as the incidence or yaw is altered. The 
method of lines has been used to provide resolution of the thin entropy layer in the 
neighborhood of the conical body and a study was also made in an attempt to provide ver- 
ification for the lift-off of the nodal singularity. Computations were made up to the incip- 
ient lift-off condition but could not be carried to larger angles of attack because compu- 
tational difficulties are encountered as the incipient lift-off condition is approached. 

Comparisons have been made with experiment and other computations for the com- 
pression side of conical delta wings with attached leading-edge shocks. Computations 
for wings with parabolic, circular arc, and flat -plate cross sections are in good agree- 
ment with other calculations and with experiment except in the region of the leading edge 
where the viscous interaction effects are most important. Over the rest of the wing the 
calculated pressures fall at the lower limit of the experimental data. Large regions of 
supersonic cross flow occur on the conical wings but contrary to the case of circular and 
elliptic cones, this condition provides no limitation and the method is applicable with no 
difficulties. Some cross-flow streamline computations are presented to illustrate the 
flow pattern on the windward side of conical wings. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., August 6, 1971. 
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APPENDIX A 


GEOMETRICAL RELATIONS 


In order to make use of the body-oriented conical coordinate system ?j,t, it is 
necessary to relate the arc length t (on the intersection of the unit sphere with the con- 
ical body) to the rectangular body coordinates x 0 ,y D . Moreover, the direction cosines 
between the body-oriented coordinates and the Cartesian coordinates X, Y, and Z are 
required to relate the velocity components and the coordinates in the two systems. These 
relations have previously been given in reference 43 but are presented here in the notation 
of this paper for the sake of completeness. 


Arc Length 

The element of length ds in the r),T, r coordinates and the Cartesian coordinates 
is given by 

(ds) 2 = r^|(d?7) 2 + (hdr)^ + (dr) 2 = (dX) 2 + (dY) 2 + (dZ) 2 


where rh is the scale factor of the r coordinate. On the conical body the rectangular 
coordinates x 0 ,y 0 are related to the Cartesian coordinates by 


*o = 

y 0 = < Y / z )o 


(Al) 


where the subscript o denotes values on the conical body. In the body-oriented coor- 
dinates, r] = 0 on the body and h = 1 at r) = 0 since t is the arc length on the con- 
ical body at r = 1. Thus, from the expression for the element of length 


dr = ± 




(A2) 


In order to integrate equation (A2), it is necessary to find relations between (— ) , 

\dx/ o 

, dX Q and the differentials dx Q and dy Q in terms of the geometric parameters 
o 

which describe the conical body. Differentiation of equations (Al) gives 


dZ 

dX 
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APPENDIX A - Continued 



The intersection of the spherical surface r — 1 with the body, in Cartesian coordinates, 
is 


X 2 + Y 2 
A o + *o 


+ Z„ - 1 = 0 


(A4a) 


or 


Z -A' l/2 
L o A 1 


(A4b) 


where Aj = 1 + x Q 2 + y o 2 . Differentiation of equation (A4) gives 

■.*> 41 ). *(§)." 

Differentiation of the equation defining the conical body, G(x 0 ,y Q ) = 0, together with 
equations (A3) gives 


(A 5) 



+ y G 
y o y 



= 0 


(A6) 


where the subscripts on G denote derivatives with respect to the indicated argument. 
Substitution from equations (A3), (A4), (A5), and (A6) into equation (A2) gives 


dr = -A^A^G" 1 dx 0 = A-U^G' 1 dy Q 


(A7a) 


d| = k T dr 


(A7b) 


where 


-^2 G x + G y + ( x o G x + ^o G y) 

Integration of equations (A7) provides the relation between the body rectangular coordi- 
nates x Q ,y 0 , the arc length t, and the transformed coordinate The integration is 
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APPENDIX A - Continued 


first carried out to determine t(N) and £(N). The line spacing is then specified in 
terms of £ and subsequent integration of equations (A7) establishes the values of x 0 
and y 0 at each line. 


Direction Cosines 

The unit vector along a ray of the body is 


®r,o = A 1 2 (**o + iy 0 + 6 ) 


and the unit normal to the body e is 

V > o 


e _ VG _ a -1/2 

7? >° |vg| 2 


"iG x + jG y - k(x 0 G x + y 0 G y ) 


where 


V = 


9X 0 dY Q az 0 


Then the unit vector e T normal to the plane r = Constant is 


e T = 


e x e 
r,o 77,0 


The unit vectors e r and e^ are found by a rotation of e r Q and e^ Q through the 
angle 77 in the plane t = Constant. Then the set of unit vectors e , § T , and § r is 
related to the Cartesian unit vectors by 


s 


a ll 

a 12 

a 13 


= 

a 21 

a 22 

a 23 

fr_ 


a 31 

a 32 

a 33 


i 

I 

k 


where 



APPENDIX A - Concluded 


a H = -x Q A^ 1 / 2 sin 7 7 + Ag^Gjj cos 77 
a 12 = -y 0 A i 1 ^ 2 sin 77 + Ag^^Gy cos 77 
a 13 = -AJ*/ 2 sin 7 ] - Ag^Ag cos 77 


a 21 " ~ ( A 1 A 2 ) ^ 2 ( G y + y 0 A 3 ) 
a 22 - ( A 1 A 2 )' 1/2 (°x +x o A 3 ) 
a 23 - ( A l A 2 )' 1 / 2 ( x o G y - y 0 G x) 


a 31 = x q A^/ 2 cos 77 + Ag^‘ //2 G X sin 77 
a 32 = y 0 A i 1//2 cos 77 + Ag 1 / 2 Gy sin 77 
a 33 = A i 1/2 cos 77 - Ag^Ag sin 77 
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APPENDIX B 


FLOW PROPERTIES BEHIND SHOCK WAVE 

Once the shock wave shape is prescribed the flow properties behind the shock can 
be evaluated from the shock jump conditions. In the case of a wing with an attached 
leading-edge shock, the flow is completely determined at the leading edge from these 
relations since the shock shape is known from the geometry there and from the stream 
conditions. 


Shock Conditions 

A set of unit vectors is constructed normal and tangential to the shock to obtain 
relations for the pressure, density, and velocity components behind the shock. Let 
n^ be the unit vector along a ray at the shock surface, n 2 be the outward unit normal 
to the shock, rig = n^ X fig, and a be the angle between the shock normal and the unit 



Figure 38.- Geometric relations between unit vectors 
at shock in a plane normal to a conical ray. 

vector e^ (fig. 38). The angle ct is given by 


tan a = 


l^s 

h dr 


h d£ 

where the scale factor h is evaluated at ?) = r) Q . This set of unit vectors is related to 
the set e T ,e ,§ r by 

hi = e r 

n 2 = cos a + e T sin a 
hg = -e sin a + § T cos a 
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APPENDIX B - Continued 


The stream velocity vector V M is taken of unit magnitude here for convenience. 
This vector, expressed in terms of the three sets of unit vectors, is 

Voo = j sin a + k cos a 


V r e r + + V T e T 


= Vi7?i+V 2 j?2 + V 3^3 


where 


V r = ^32 s ^ n a + a 33 cos a 


~ a 12 


" a 22 




sin a + a ^2 
sin a + ag 3 


cos a 
cos a 


cos a + V T sin a 
sin <j + V T cos a 


and a is the angle of incidence. The direction cosines a^j between the vectors 
e^,§ T ,e r and I,j,k are given in appendix A. The stream velocity components 
and Vg are tangential to the shock and Vg is normal to the shock. The shock 
angle /3 then is given by 

sin /3 = -n 2 • V M = -V 2 

and the pressure and density for an ideal gas are given by (ref. 59) 

p = j2yM oo 2 sin 2 /3 - (y - 1)1 

y(y + ^M^c* 2 

(y + l)M 00 2 sin% 

p = 

2 + (y - ^M^sin 2 ^ 
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APPENDIX B - Continued 


From the Rankine-Hugoniot relations, the required velocity components behind the shock, 
referenced to the critical speed, are 



v = -Jjr(p -1 V2 cos a " V 3 sin ° 
w = sin cr + V 3 cos cr) 


Attached Shock at Wing Leading Edge 

For the computation of the flow about wings with an attached leading -edge shock, 
the angle a at the leading edge can be computed from the stream conditions and the 
wing geometry. Once a is evaluated, all the flow quantities at the leading edge are 
computed directly from the foregoing shock relations. The wing leading edge lies in the 
X,Z -plane and makes an angle 9 with the Z-axis ^the sweep angle A = ^ - 0^. To 

determine a, it is only necessary to consider the flow quantities in a plane normal to the 
wing leading edge as shown in figure 39. The subscript n identifies values in this 
plane . 


Su rface 

Figure 39 -- Geometric relations at wing leading edge in a plane normal to edge. 

The component of the stream velocity vector (v m is of unit magnitude) in the plane 
normal to the wing leading edge, V n , and its magnitude, V n , are 

V n = -I sin 6 cos 9 cos a + j sin a + k sin 2 0 cos a 

V n = (sin 2 0 cos 2 or + sin^a;)*'^ 



The free -stream Mach number in this plane is M n = M 00 V n< 

Let j u be the opening angle of the wing leading edge (measured from the X-axis) 
in a plane Z = Constant. Then p is computed from the geometry at the leading edge 
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APPENDIX B — Concluded 


by 


tan (J. = 


d y 0 

dx„ 


and the angle p n is found from 

, tan u 

tan u n = £ 

n cos 9 


From the wing geometry, the vector e^ is 

e^ = I sin ju n cos 9 - j cos ju n - k sin sin 9 
and 5 n is given by 


sin 


V e 1 

5 n = - — -L = V n ' A (sin 9 cos a sin p n + sin a cos |U n ^ 


The shock angle /3 n is found from the equation relating the flow deflection 5 n and /3 n 
(ref. 59) 


cot 6 n = tan 0 n 


(r + DM n 2 


2(M n 2 sin 2 )3 n - l) 


and a = - S n . The flow quantities at the wing leading edge are computed from the 

shock relations with this value of a and replaced by M n . 
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APPENDIX C 


FORCE AND MOMENT COEFFICIENTS 

Once the surface values of the pressure have been evaluated, the force and moment 
coefficients can be computed by numerical integration. The unit vector normal to the 
surface e (appendix A) is 

l/ 9 U 



and the element of area on the conical surface is r dr dr. The normal force F n on 
this elemental area is 

dF\ o 

2— = -e c* ( p - p„A r dr dr 

o y 2 ' °°) 

‘oo v oo 

= -%0 C+2 (P -Poo)^ 1 r d ? dr 

where p^ = ^ Then the components of the force along the Y and Z direc- 

tions are 

d.F 

Y_ = _c*2j p _ Poo ) A -l /2 Gy |;l r d£ dr (Cl) 

Poo^oo 

= c*2(p . p„) A 2 1/2 (x 0 G x + y^S^r d| dr (C2) 

Poo' oo 

The integration is carried out over the conical body to the surface Z D = 1. Now 

r = ZqAj* 7 ; thus, at Z Q = 1, r = . The force coefficients then, after integration 

on r, are 

C Y = ^f(»-^) A l A i 1/2 V> 
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APPENDIX C - Continued 


£ 

C Z - ^ U(P ' p «) a i ^ 1/2 ( x ^ * y » G y) { / d| 

where the upper limit £ u = | N for elliptic cones and £ u = £ N+1 for wing configurations 
and A is the reference area. The drag and lift coefficients are determined from 

Cd = Cy sin a. + cos a 


Cl = Cy cos a - C^ sin a 


Expressions for the differential moments, or the center -of -pressure location mea- 
sured from the apex, are found by multiplying equation (Cl) by Z 0 and equation (C2) 
by Y 0 . Thus, with 


z o = rA l 


1/2 


Y = v Z 
o y o o 


the coordinates of the center of pressure Y C p and Z^, after integration on r, are 


Y„„ = 


4c*^ 


cp 3AC Z 


r = 4c*^ 

J CP " 3ACt 


^ “(P - P-) A l A 2 1/2y o( X o G X + y 0 G y) f r ld « 

^ U (P-P») A 1 A 2 I/ V>=I 


and the moment about the X-axis is 


C m,x = Y cp C Z - Z cp C Y 


For the elliptic cone the reference area A is taken as the base area 
A = 7r^)tan 2 e 

and for the wing the reference area is the plan area 


A = tan 9 
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APPENDIX C - Concluded 


The integrals for the force coefficients and the coordinates of the center of pressure can 
be evaluated with the use of Simpson’s rule if the number of increments At is even; 
hence, it is required that the number of lines N is odd for elliptic cones and even for 
the wings with attached leading-edge shocks. 


1 3 



REFERENCES 


1. Busemann, A.: Aerodynamic Lift at Supersonic Speeds. Ae. Techl. 1201, Rep. 

No. 2844, Brit. A.R.C., Feb. 3, 1937. (From Luftfahrtforschung, Bd. 12, Nr. 6, 
Oct. 3, 1935, pp. 210-220. 

2. Jones, D. J.: Numerical Solutions of the Flow Field for Conical Bodies in a Super- 

sonic Stream. Aeronaut. Rep. LR-507 (NRC No. 10361), Nat. Res. Counc. Can. 
(Ottawa), July 1968. 

3. Ndefo, D. Ejike: A Numerical Method for Calculating Steady Unsymmetrical Super- 

sonic Flow Past Cones. Rep. No. AS-69-11 (AFOSR Grant 268-68), U.S. Air 
Force, May 1969, (Available from DDC as AD 691 270.) 

4. South, Jerry C., Jr,; and Klunker, E. B. : Methods for Calculating Nonlinear Conical 

Flows. Analytic Methods in Aircraft Aerodynamics, NASA SP-228, 1970, 
pp. 131-155; Discussion, pp. 156-158. 

5. Stone, A. H. : On Supersonic Flow Past a Slightly Yawing Cone. J. Math. Phys., 

vol. XXVH, No. 1, Apr. 1948, pp. 67-81. 

6. Staff of Comput. Section, Center of Anal. (Under dir. of Zdenek Kopal): Tables of 

Supersonic Flow Around Yawing Cones. Tech. Rep. No. 3 (NOrd Contract 
No. 9169), Massachusetts Inst. Technol., 1947. 

7. Taylor, G. I.; and Maccoll, J. W.: The Air Pressure on a Cone Moving at High 

Speeds. Proc. Roy. Soc. (London), ser. A, vol. 139, no. 838, Feb. 1, 1933, 
pp. 278-311. 

8. Ferri, Antonio: Supersonic Flow Around Circular Cones at Angles of Attack. NACA 

Rep. 1045, 1951. (Supersedes NACA TN 2236.) 

9. Melnik, R. E.: Vortical Layers in Supersonic Conical Flow. Hypersonic Boundary 

Layers and Flow Fields. AGARD CP No. 30, 1968, pp. 26-1 - 26-20. 

10. Woods, B. A.: The Supersonic Flow Past an Elliptic Cone. Aeronaut. Quart., 

vol. XX, pt. 4, Nov. 1969, pp. 382-404. 

11. Munsen, A. G.: A Uniformly Valid Solution for the Flow Over an Inclined Cone Using 

the Method of Matched Asymptotic Expansions. AFOSR 64-1056, U.S. Air Force, 
Apr. 1964. (Available from DDC as AD 603 154.) 

12. Moretti, Gino: Inviscid Flowfield About a Pointed Cone at an Angle of Attack. AIAA 

J., vol. 5, no. 4, Apr. 1967, pp. 789-791. 

13. Babenko, K. I.; Voskresenskii, G. P.; Lyubimov, A. N.; and Rusanov, V. V.: Three- 

Dimensional Flow of Ideal Gas Past Smooth Bodies. NASA TT F-380, 1966. 


74 



14. Rakich, John V.: A Method of Characteristics for Steady Three-Dimensional Super- 

sonic Flow With Application to Inclined Bodies of Revolution. NASA TN D-5341, 
1969. 

15. Babenko, K. I.; and Rusanov, V. V.: Difference Methods of Solving Three- 

Dimensional Problems in Gas Dynamics. NASA TT F-10,827, 1967. 

16. Gonidou, Rene: Supersonic Flows Around Cones at Incidence. NASA TT F-11,473, 

1968. 

17. Voskresenskii, G. P.: Chislennoe Reshenie Zadachi Obtekaniia Proizvol’noi 

Poverkhnosti Treugul'nogo Kryla v Oblasti Szhatiia Sverkhzvukovym Potokom 
Gaza (Numerical Solution of the Problem of a Supersonic Gas Flow Past an 
Arbitrary Surface of a Delta Wing in the Compression Region). Izv. Akad. Nauk 
SSSR, Mekh. Zhidk. Gaza, no. 4, 1968, pp. 134-142. 

18. Powers, S. A.; and Beeman, E. R., Jr.: Flow Fields Over Sharp Edged Delta Wings 

With Attached Shocks. NASA CR-1738, 1971. 

19. Kutler, Paul: Application of Selected Finite Difference Techniques to the Solution 

of Conical Flow Problems. Ph. D. Diss., Iowa State Univ., 1969. 

20. Kutler, Paul; and Lomax, Harvard: A Systematic Development of the Supersonic 

Flow Fields Over and Behind Wings and Wing-Body Configurations Using a Shock- 
Capturing Finite-Difference Approach. AIAA Paper No. 71-79, Jan 1971. 

21. Busemann, A.: Driicke auf Kegelformige Spitzen bei Bewegung mit Uberschall- 

geschwindigkeit. Z. Angew. Math. Mech., Bd. 9, Heft 6, Dec. 1929, pp. 496-498. 

22. Hayes, Wallace D.; and Probstein, Ronald F.: Hypersonic Flow Theory. Vol. I — 

Inviscid Flows. Second ed., Academic Press, 1966. 

23. Briggs, Benjamin R.: The Numerical Calculation of Flow Past Conical Bodies Sup- 

porting Elliptic Conical Shock Waves at Finite Angles of Incidence. NASA 
TN D-340, 1960. 

24. Mauger, F. E.: Steady Supersonic Flow Past Conical Bodies. A.R.D.E. Rep. (B)3/60, 

Brit. War Office, May 1960. 

25. Stocker, P. M.; and Mauger, F. E.: Supersonic Flow Past Cones of General Cross- 

Section. J. Fluid Mech., vol. 13, pt. 3, July 1962, pp. 383-399. 

26. Eastman, E. W.; and Omar, M. E.: Flow Fields About Highly Yawed Cones by the 

Inverse Method. AIAA J., vol. 3, no. 9, Sept. 1965, pp. 1782-1784. 

27. Makhin, N. A.; and Syagayev, V. F.: Numerical Solution of the Problem of Supersonic 

Flow Past Conical Bodies at an Angle of Attack. Fluid Dyn., vol. 1, no. 1, Jan. -Feb. 
1966, pp. 100-101. 


75 



28. Maslen, Stephen H.: Supersonic Conical Flow. NACA TN 2651, 1952. 

29. Babayev, D. A.: Numerical Solution of the Problem of Flow Round the Upper Surface 

of a Triangular Wing by a Supersonic Stream. U.S.S.R. Comput. Math. Math. Phys., 
no. 2, 1963, pp. 296-308. 

30. Babaev, D. A.: Numerical Solution of the Problem of Supersonic Flow Past the 

Lower Surface of a Delta Wing. AIAA J., vol. 1, no. 9, Sept. 1963, pp. 2224-2231. 

31. Dorodnitsyn, A. A.: A Contribution to the Solution of Mixed Problems of Transonic 

Aerodynamics. Vol. 2 of Advances in Aeronautical Sciences, Pergamon Press, 

Inc., 1959, pp. 832-844. 

32. Belotserkovskiy, O. M.: Supersonic Gas Flow Around Blunt Bodies - Theoretical 

and Experimental Investigations. NASA TT F-453, 1967. 

33. South, Jerry C., Jr.: Application of the Method of Integral Relations to Supersonic 

Nonequilibrium Flow Past Wedges and Cones. NASA TR R-205, 1964. 

34. Chushkin, P. I.; and Shchennikov, V. V. (B. A. Woods, transl.): Calculation of Certain 

Conical Flows Without Axial Symmetry. Lib. Transl. No. 926, Brit. R.A.E., Dec. 
1960. 

35. Brook, John W.: The Method of Integral Relations for Conical Flow - Theoretical 

Analysis. Res. Mem. RM-193 (Contract No. AF 33(616) -6400), Grumman Aircraft 
Eng. Corp., Oct. 1961. (Available from DDC as AD 266 501.) 

36. Liskovets, O. A.: The Method of Lines (Review). Differential Equations, vol. 1, 

no. 12, Dec. 1965, pp. 1308-1323. 

37. Roslyakov, G. S.; and Telenin, G. F.: Survey of Computational Work on Steady 

Axisymmetric Gas Flow Carried out at the Computational Center of Moscow State 
University. Numerical Methods in Gas Dynamics, G. S. Roslyakov and L. A. 

Chudov, eds., NASA TT F-360, Israel Program Sci. Transl., 1966, pp. 1-11. 
(Available from CFSTI, U.S. Dep. Com.) 

38. Telenin, G. F.; and Tinyakov, G. P.: A Method of Calculating the Three-Dimensional 

Flow Past a Body With an Attached Shock Wave. Sov. Phys. - Doklady, vol. 9, 
no. 2, Aug. 1964, pp. 132-133. 

39. Tiniakov, G. P.; Investigation of the Three-Dimensional Supersonic Flow Around 

Ellipsoids of Revolution. Lockheed Missiles and Space Co. Transl. (From Izv. 
Akad. Nauk SSSR, Otd. Tekhn, Nauk, Mekhan, i Mashinostr., no. 6, 1965, 
pp. 10-19.) 

40. Makhin, N. A.; and Syagayev, V. F.: Numerical Solution of the Problem of Supersonic 

Flow at an Angle of Attack Past Conical Bodies. NASA TT F-10,481, 1966. 


76 



41. Bazzhin, A. P.; and Chelysheva, I. F.: Primenenie Metoda Priamykh k Raschetu 

Obtekaniia Konicheskikh Tel pri Bol'shikh Uglakh Ataki (Application of the Straight- 
Line Method in the Calculation of Flows Past Conical Bodies at Large Angles of 
Attack]. Izv. Akad. Nauk SSSR, Mekh. Zhidk. Gaza, no. 3, 1967, pp. 119-123. 

42. Bazzhin, A. P.; Trusova, O. N.; and Chelysheva, I. F.: Raschet Techenii 

Sovershennogo Gaza Okolo Ellipticheskikh Konusov pri Bol'shikh Uglakh Ataki 
[Calculation of Ideal -Gas Flows Around Elliptical Cones at Large Angles of Attack], 
Izv. Akad. Nauk SSSR, Mekh. Zhidk. Gaza, no. 4, 1968, pp. 45-51. 

43. Scheuing, Richard A.; Brook, John W.; Mead, Harold R.; Melnik, Robert E.; Hayes, 

Wallace D.; Donaldson, Coleman duP.; Gray, K. Evan; and Sullivan, Roger D.: 
Theoretical Prediction of Pressures in Hypersonic Flow With Special Reference to 
Configurations Having Attached Leading-Edge Shock - Pt. I. Theoretical 
Investigation. ASD TR 61-60, Pt. I, U.S. Air Force, May 1962. 

44. McConnell, A. J.: Applications of Tensor Analysis. Dover Publ., Inc., 1957. 

45. Nielsen, Kaj L.: Methods in Numerical Analysis. Macmillan Co., c.1956. 

46. Mysovskih, I. P. (L. B. Rail, trans.): Lectures on Numerical Methods. Walters- 

Noordhoff Publ. (Groningen, Netherlands), c.1969. 

47. Hord, Richard A.: An Approximate Solution for Axially Symmetric Flow Over a Cone 

With an Attached Shock Wave. NACA TN 3485, 1955. 

48. Hadamard, Jacques: Lectures on Cauchy's Problem in Linear Partial Differential 

Equations. Dover Publ., Inc., c.1952. 

49. Garabedian, P. R.: Partial Differential Equations. John Wiley & Sons, Inc., c.1964. 

50. Jones, D. J.: Tables of Inviscid Supersonic Flow About Circular Cones at Incidence 

y = 1.4. AGARDograph 137, Pts. I and II, Nov. 1969. 

51. Tracy, Richard R.: Hypersonic Flow over a Yawed Circular Cone. Hypersonic Res. 

Proj. Mem. No. 69 (Contract No. DA-31- 124 -ARO(D) -33), Graduate Aeronaut. Lab., 
California Inst. Technol., Aug. 1, 1963. 

52. Melnik, R. E.: Vortical Singularities in Conical Flow. AIAA J., vol. 5, no. 4, Apr. 

1967, pp. 631-637. 

53. Babenko, K. I.: Investigation of a Three-Dimensional Supersonic Gas Flow Around 

Conic Bodies. Applied Mechanics, Henry Gortler and Peter Sorger, eds.. 

Springer -Verlag, 1966, pp. 749-755. 

54. Chapkis, Robert L. : Hypersonic Flow Over an Elliptic Cone: Theory and Experiment. 

J. Aerosp. Sci., vol. 28, no. 11, Nov. 1961, pp. 844-854. 


77 



55. Martellucci, Anthony: An Extension of the Linearized Characteristics Method for 

Calculating the Supersonic Flow Around Elliptic Cones. J. Aerosp. Sci., vol. 27, 
no. 9, Sept. 1960, pp. 667-674. 

56. Chiang, C. W.; and Wagner, Richard D., Jr.: Analysis of Supersonic Conical Flows. 

NASA TN D-5884, 1970. 

57. Vincenti, Walter G.; and Fisher, Newman H., Jr.: Calculation of the Supersonic 

Pressure Distribution on a Single -Curved Tapered Wing in Regions Not Influenced 
by the Root or Tip. NACA TN 3499, 1955. 

58. Mead, Harold R.; and Koch, Frank: Theoretical Prediction of Pressures in Hyper- 

sonic Flow With Special Reference to Configurations Having Attached Leading- 
Edge Shock - Pt. n. Experimental Pressure Measurements at Mach 5 and 8. 

ASD TR 61-60, Pt. II, U.S. Air Force, May 1962. 

59. Ames Research Staff: Equations, Tables, and Charts for Compressible Flow. NACA 

Rep. 1135, 1953. (Supersedes NACA TN 1428.) 


78 


NASA-Langley, 1971 12 L-7813 


NATIONAL AERONAUTICS AND SPACE ADMISTRATION 
WASHINGTON, D.C. 20546 

official business FIRST CLASS MAIL 

penalty for private USE $300 


POSTAGE AND FEES PAID 
NATIONAL AERONAUTICS AND 
SPACE ADMINISTRATION 



020 001 Cl U 12 711001 S00903DS 
DEPT OF THE AIB FORCE 
AF SYSTEMS COMMAND 
AF WEAPONS LAB (WLOL) 

ATTN: E LOU BOWMAN, CHIEF TECH LIBRARY 
KIRTLAND AFB NM 87117 


a ctcd . ^ Undeliverable ( Section 158 

POSTMASTER. postal Manual ) Do Noc Return 


* 


"The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of hutnan knowl- 
I - edge of phenomena in the atynosphere and space. The Administration 

shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof .” 

— National Aeronautics and Space Act of 1958 

NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. . 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs 0 
Technology Utilization Reports and 
Technology Surveys. 


s Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION OFFICE 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 



I 

V 

i 


I 


/4 

1 


a 



